text:
from __future__ import unicode_literals
from __future__ import absolute_import

# 3rd party modules
import numpy as np
from datetime import datetime
import webbrowser

# QGIS modules
from qgis.core import *
from qgis.PyQt.QtCore import Qt
from qgis.PyQt.QtWidgets import *
from qgis.PyQt import uic

# promaides modules
from .interpolate import RasterInterpolator
from .environment import get_ui_path
from .version import *
from .utils import *

try:
    from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as Canvas
    from matplotlib.figure import Figure
    from matplotlib.widgets import Cursor
    PREVIEW_ENABLED = True
except ImportError:
    PREVIEW_ENABLED = False
    Canvas, Figure, Cursor = None, None, None

UI_PATH = get_ui_path('ui_river_profile_export.ui')

# This plugin exports the 1D-river file for the HYD-module of ProMaIdes from a line shape file;
# A name field (string) is required within the line layer
class PluginDialog(QDialog):

    def __init__(self, iface, parent=None, flags=Qt.WindowFlags()):
        QDialog.__init__(self, parent, flags)
        uic.loadUi(UI_PATH, self)

        self.iface = iface
        self.input_layer = None
        self.raster_layer = None
        self.roughness_layer = None
        #add channel info
        self.channel_info = None

        # self.name_box.setExpression('')
        # self.station_box.setExpression('')
        self.delta_box.setExpression('0.1')
        self.connection_box.setExpression("'standard'")
        self.type_box.setExpression("'river'")
        self.initial_box.setExpression('0.5')
        self.point_bc_enabled_box.setExpression('False')
        self.point_bc_stationary_box.setExpression('True')
        self.point_bc_value_box.setExpression('0')
        self.lateral_bc_enabled_box.setExpression('False')
        self.lateral_bc_stationary_box.setExpression('True')
        self.lateral_bc_value_box.setExpression('0')
        self.overflow_left_enabled_box.setExpression('True')
        self.overflow_left_poleni_box.setExpression('0.577')
        self.overflow_right_enabled_box.setExpression('True')
        self.overflow_right_poleni_box.setExpression('0.577')
        self.localbridgeheight_box.setExpression('1.5')
        self.bridgebodyheight_box.setExpression('4.0')

        self.groupBox_bridge.setCollapsed(True)
        self.mGroupBox_6.setCollapsed(True)
        self.mGroupBox.setCollapsed(True)
        self.mGroupBox_4.setCollapsed(True)
        self.mGroupBox_3.setCollapsed(True)
        self.mGroupBox_2.setCollapsed(True)
        self.groupBox_main.setCollapsed(True)
        self.mGroupBox_4.setCollapsed(False)

        #set DGM layer
        self.raster_layer_box.setFilters(QgsMapLayerProxyModel.RasterLayer)
        self.raster_layer_box.layerChanged.connect(self.setRasterLayer)
        self.raster_band_box.setEnabled(False)

        self.method_box.addItem('nearest neighbor (downscaling/upscaling)')
        self.method_box.addItem('bi-linear (downscaling)')
        self.method_box.addItem('bi-cubic (downscaling)')

        #set roughness layer
        self.roughness_layer_box.setFilters(QgsMapLayerProxyModel.RasterLayer)
        self.roughness_layer_box.layerChanged.connect(self.setRoughnessLayer)
        self.roughness_layer_box.setLayer(None)

        #set channel layer
        self.comboBox_main.setFilters(QgsMapLayerProxyModel.PolygonLayer)
        self.comboBox_main.layerChanged.connect(self.setChannelLayer)
        self.comboBox_main.setLayer(None)

        self.browse_button.clicked.connect(self.onBrowseButtonClicked)
        self.HelpButton.clicked.connect(self.Help)
        self.iface.currentLayerChanged.connect(self.setInputLayer)

        #############################################
        #autostationing button function
        def autoclicked(state):
            if state > 0:
                self.profileid_box.setEnabled(True)
                self.station_box.setEnabled(False)
                self.preview_button.setDisabled(True)
            else:
                self.profileid_box.setEnabled(False)
                self.station_box.setEnabled(True)
                self.preview_button.setDisabled(False)

        self.profileid_box.setEnabled(False)
        self.autostation_box.stateChanged.connect(autoclicked)
        #############################################

        if PREVIEW_ENABLED:
            self.figure = Figure(figsize=(10, 4), dpi=100)
            self.axes = self.figure.add_subplot(111)
            self.axes.set_xlabel('Stationing')
            self.axes.set_ylabel('Elevation [m+NHN]')
            self.figure.tight_layout()
            self.canvas = Canvas(self.figure)
            self.canvas.setParent(self)
            self.canvas.setSizePolicy(QSizePolicy.Expanding, QSizePolicy.Expanding)
            self.layout().addWidget(self.canvas)
            self.canvas.setVisible(False)
            self.cursor = None
            self.preview_button.toggled.connect(self.canvas.setVisible)
            self.preview_button.toggled.connect(self.update_button.setEnabled)
            self.update_button.clicked.connect(self.updateLongitudinalPreview)
        else:
            self.preview_button.setDisabled(True)
            self.figure = None
            self.axes = None
            self.canvas = None
            self.cursor = None

        self.setInputLayer(self.iface.activeLayer())
        self.setRasterLayer(self.raster_layer_box.currentLayer())

    def __del__(self):
        #self.iface.currentLayerChanged.disconnect(self.setInputLayer)
        print('Done')

    def Help(self):
        webbrowser.open("https://promaides.myjetbrains.com/youtrack/articles/PMDP-A-45/1D-River-Profile-Export")

    def onBrowseButtonClicked(self):
        current_filename = self.filename_edit.text()
        new_filename, __ = QFileDialog.getSaveFileName(self.iface.mainWindow(), '1D-River Profile Export', current_filename)
        if new_filename != '':
            self.filename_edit.setText(new_filename)
            self.filename_edit.editingFinished.emit()

    def setInputLayer(self, layer):
        """
        """
        if not layer:
            self.input_layer = None
            self.input_label.setText(u'No layer selected.
'
                                     u'Please select a polyline layer.')
        else:
            layer_name = layer.name()
            if layer.type() == QgsMapLayer.VectorLayer:
                if layer.geometryType() == QgsWkbTypes.LineGeometry:
                    self.input_layer = layer
                    if layer.selectedFeatureCount():
                        self.input_label.setText(u'Input layer is "{}" with {} selected feature(s).'
                                                 .format(layer_name, layer.selectedFeatureCount()))
                    else:
                        self.input_label.setText(u'Input layer is "{}" with {} feature(s).'
                                                 .format(layer_name, layer.featureCount()))

                else:
                    self.input_layer = None
                    self.input_label.setText(u'Selected layer "{}" is not a polyline layer.
'
                                             u'Please select a polyline layer.'
                                             .format(layer_name))
            else:
                self.input_layer = None
                self.input_label.setText(u'Selected layer "{}" is not a vector layer.
'
                                         u'Please select a polyline layer.'
                                         .format(layer_name))

        self.name_box.setLayer(self.input_layer)
        self.station_box.setLayer(self.input_layer)
        self.delta_box.setLayer(self.input_layer)
        self.connection_box.setLayer(self.input_layer)
        self.type_box.setLayer(self.input_layer)
        self.initial_box.setLayer(self.input_layer)
        self.profileid_box.setLayer(self.input_layer)
        self.point_bc_enabled_box.setLayer(self.input_layer)
        self.point_bc_stationary_box.setLayer(self.input_layer)
        self.point_bc_value_box.setLayer(self.input_layer)
        self.lateral_bc_enabled_box.setLayer(self.input_layer)
        self.lateral_bc_stationary_box.setLayer(self.input_layer)
        self.lateral_bc_value_box.setLayer(self.input_layer)
        self.overflow_left_enabled_box.setLayer(self.input_layer)
        self.overflow_left_poleni_box.setLayer(self.input_layer)
        self.overflow_right_enabled_box.setLayer(self.input_layer)
        self.overflow_right_poleni_box.setLayer(self.input_layer)
        self.localbridgeheight_box.setLayer(self.input_layer)
        self.bridgebodyheight_box.setLayer(self.input_layer)

        self.updateButtonBox()

    def setRasterLayer(self, layer):
        self.raster_layer = layer
        if layer:
            self.raster_band_box.setEnabled(True)
            self.raster_band_box.setRange(1, layer.bandCount())
        else:
            self.raster_band_box.setEnabled(False)

        self.updateButtonBox()

    def setRoughnessLayer(self, layer):
        self.roughness_layer = layer
        if layer:
            self.roughness_band_box.setEnabled(True)
            self.roughness_band_box.setRange(1, layer.bandCount())
        else:
            self.roughness_band_box.setEnabled(False)

    def setChannelLayer(self, layer):
            self.channel_info = layer

    def updateButtonBox(self):
        if self.input_layer and self.raster_layer:
            self.button_box.button(QDialogButtonBox.Ok).setEnabled(True)
        else:
            self.button_box.button(QDialogButtonBox.Ok).setEnabled(False)

    def updateLongitudinalPreview(self):
        if not PREVIEW_ENABLED:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Preview is disabled since Python module "matplotlib" is not installed!'
            )
            return

        flip_directions = self.flip_directions_box.isChecked()
        addfullriver = self.fullriver_box.isChecked()
        autostation = self.autostation_box.isChecked()
        abs_init = self.abs_init_box.isChecked()
        adjust_elevation = self.adjust_elevation_box.isChecked()

        input_layer = self.input_layer
        if not input_layer:
            return

        if input_layer.selectedFeatureCount():
            selected = True
            features = input_layer.selectedFeatures()
        else:
            selected = False
            features = input_layer.getFeatures()

        stations, ok = QgsVectorLayerUtils.getValues(input_layer, self.station_box.expression(), selected)
        if not autostation:
            if not ok:
                self.iface.messageBar().pushCritical(
                    'River Profile Export',
                    'Invalid expression for stations!'
                )
                return

        init_values, ok = QgsVectorLayerUtils.getValues(input_layer, self.initial_box.expression(), selected)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for initial condition!'
            )
            return

        dem_layer = self.raster_layer
        dem_band = self.raster_band_box.value()
        dem_method = self.method_box.currentText()
        dem_nan = self.nan_box.value()
        dem_interpol = RasterInterpolator(dem_layer, dem_band, 10, 10, dem_method, dem_nan)

        base, left, right, h = [], [], [], []

        # iterate over profiles and extract base elevation,
        # left and right bank elevation, and init value
        for i, f in enumerate(features):

           # line = f.geometry().asPolyline()
            line = []
            buff = f.geometry()
            for p in buff.vertices():
                line.append(p)

            if flip_directions:
                line = list(reversed(line))

            z = [dem_interpol(QgsPointXY(p)) for p in line]

            if adjust_elevation:
                if z[0] <= z[1]:
                    z[0] = z[1] + 0.01
                if z[-1] <= z[-2]:
                    z[-1] = z[-2] + 0.01

            i_min = np.argmin(z)
            # minimum profile elevation
            z_min = z[int(i_min)]

            base.append(z_min)
            left.append(z[0])
            right.append(z[-1])

            init_value = float(init_values[i])
            # init values are absolute values
            if not abs_init:
                init_value += z_min

            h.append(init_value)

        s, base, left, right, h = list(zip(*sorted(zip(stations, base, left, right, h))))

        self.axes.cla()
        self.axes.plot(s, base, 'k', lw=1.5, label='base')
        self.axes.plot(s, left, 'r', lw=1.0, label='left bank')
        self.axes.plot(s, right, 'g', lw=1.0, label='right bank')
        self.axes.plot(s, h, 'b--', lw=0.5, label='init. h')
        self.axes.legend()
        self.axes.set_xticks(s)
        self.axes.set_xticklabels(list(map(str, s)), rotation=90, fontdict=dict(fontsize=9))
        self.axes.set_xlabel('Stationing')
        self.axes.set_ylabel('Elevation [m+NHN]')
        self.figure.tight_layout()
        self.cursor = Cursor(self.axes, useblit=True, color='0.5', linewidth=0.5)
        self.canvas.draw()

class RiverProfileExport(object):

    def __init__(self, iface):
        self.iface = iface
        self.act = None
        self.dialog = None
        self.cancel = False
        self.act = QAction('1D-River Profile Export', iface.mainWindow())
        self.act.triggered.connect(self.execDialog)

    def initGui(self, menu=None):
        if menu is not None:
            menu.addAction(self.act)
        else:
            self.iface.addToolBarIcon(self.act)

    def unload(self, menu=None):
        if menu is None:
            menu.removeAction(self.act)
        else:
            self.iface.removeToolBarIcon(self.act)

    def execDialog(self):
        """
        """
        self.dialog = PluginDialog(self.iface, self.iface.mainWindow())
        self.dialog.accepted.connect(self.execTool)
        self.dialog.rejected.connect(self.quitDialog)
        self.dialog.setModal(False)
        self.act.setEnabled(False)
        # add a filter to the combo box of the filed selection;
        self.dialog.name_box.setFilters(QgsFieldProxyModel.String)
        self.dialog.type_box.setFilters(QgsFieldProxyModel.String)
        self.dialog.connection_box.setFilters(QgsFieldProxyModel.String)
        self.dialog.point_bc_enabled_box.setFilters(QgsFieldProxyModel.String)
        self.dialog.point_bc_stationary_box.setFilters(QgsFieldProxyModel.String)
        self.dialog.lateral_bc_enabled_box.setFilters(QgsFieldProxyModel.String)
        self.dialog.lateral_bc_stationary_box.setFilters(QgsFieldProxyModel.String)
        self.dialog.overflow_right_enabled_box.setFilters(QgsFieldProxyModel.String)
        self.dialog.overflow_left_enabled_box.setFilters(QgsFieldProxyModel.String)

        self.dialog.station_box.setFilters(QgsFieldProxyModel.Numeric)
        self.dialog.delta_box.setFilters(QgsFieldProxyModel.Numeric)
        self.dialog.initial_box.setFilters(QgsFieldProxyModel.Numeric)
        self.dialog.profileid_box.setFilters(QgsFieldProxyModel.Numeric)
        self.dialog.point_bc_value_box.setFilters(QgsFieldProxyModel.Numeric)
        self.dialog.lateral_bc_value_box.setFilters(QgsFieldProxyModel.Numeric)
        self.dialog.overflow_left_poleni_box.setFilters(QgsFieldProxyModel.Numeric)
        self.dialog.overflow_right_poleni_box.setFilters(QgsFieldProxyModel.Numeric)
        self.dialog.localbridgeheight_box.setFilters(QgsFieldProxyModel.Numeric)
        self.dialog.bridgebodyheight_box.setFilters(QgsFieldProxyModel.Numeric)

        self.dialog.groupBox_bridge.setCollapsed(True)
        self.dialog.mGroupBox_6.setCollapsed(True)
        self.dialog.mGroupBox.setCollapsed(True)
        self.dialog.mGroupBox_4.setCollapsed(True)
        self.dialog.mGroupBox_3.setCollapsed(True)
        self.dialog.mGroupBox_2.setCollapsed(True)
        self.dialog.groupBox_main.setCollapsed(True)
        self.dialog.mGroupBox_4.setCollapsed(False)

        self.dialog.show()

    def scheduleAbort(self):
        self.cancel = True

    def quitDialog(self):
        self.dialog = None
        self.act.setEnabled(True)
        self.cancel = False

    def execTool(self):

        abs_init = self.dialog.abs_init_box.isChecked()
        flip_directions = self.dialog.flip_directions_box.isChecked()
        addfullriver = self.dialog.fullriver_box.isChecked()
        autostation = self.dialog.autostation_box.isChecked()
        adjust_elevation = self.dialog.adjust_elevation_box.isChecked()

        filename = self.dialog.filename_edit.text()
        if not filename:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'No output filename given!'
            )
            self.quitDialog()
            return

        input_layer = self.dialog.input_layer
        if input_layer.selectedFeatureCount():
            os = True
            features = input_layer.selectedFeatures()
            feature_count = input_layer.selectedFeatureCount()
        else:
            os = False
            features = input_layer.getFeatures()
            feature_count = input_layer.featureCount()

        names, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.name_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for names!'
            )
            self.quitDialog()
            return

        stations, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.station_box.expression(), os)
        if not autostation:
            ##################################
            #test station values to be numeric
            def isnumeric(obj):
                try:
                    obj + 0
                    return True
                except TypeError:
                    return False
            for i in stations:
                if isnumeric(i)==False:
                    self.iface.messageBar().pushCritical(
                        'River Profile Export',
                        'Invalid values for stations!'
                    )
                    self.quitDialog()
                    return
            ####################################
            if not ok:
                self.iface.messageBar().pushCritical(
                    'River Profile Export',
                    'Invalid expression for stations!'
                )
                self.quitDialog()
                return

        deltas, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.delta_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for delta table!'
            )
            self.quitDialog()
            return

        conn_types, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.connection_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for connection type!'
            )
            self.quitDialog()
            return

        profile_types, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.type_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for profile type!'
            )
            self.quitDialog()
            return

        init_values, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.initial_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for initial condition!'
            )
            self.quitDialog()
            return

        if autostation:
            profileids, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.profileid_box.expression(), os)

            #####################################
            #test profile id values to be numeric
            def isnumeric(obj):
                try:
                    obj + 0
                    return True
                except TypeError:
                    return False
            for i in profileids:
                if isnumeric(i)==False:
                    self.iface.messageBar().pushCritical(
                        'River Profile Export',
                        'Invalid values for profile ids!'
                    )
                    self.quitDialog()
                    return
            ####################################
            if not ok:
                self.iface.messageBar().pushCritical(
                    'River Profile Export',
                    'Invalid expression for profile id!'
                )
                self.quitDialog()
                return

        point_bc_flags, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.point_bc_enabled_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for point boundary enabled!'
            )
            self.quitDialog()
            return

        point_bc_stationary_flags, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.point_bc_stationary_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'Coastline Export',
                'Invalid expression for point boundary stationary!'
            )
            self.quitDialog()
            return

        point_bc_values, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.point_bc_value_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for point boundary value!'
            )
            self.quitDialog()
            return

        lateral_bc_flags, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.lateral_bc_enabled_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for lateral boundary enabled!'
            )
            self.quitDialog()
            return

        lateral_bc_stationary_flags, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.lateral_bc_stationary_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for lateral boundary stationary!'
            )
            self.quitDialog()
            return

        lateral_bc_values, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.lateral_bc_value_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for lateral boundary value!'
            )
            self.quitDialog()
            return

        overflow_left_flags, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.overflow_left_enabled_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for overflow left enabled!'
            )
            self.quitDialog()
            return

        poleni_left_values, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.overflow_left_poleni_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for left poleni factor!'
            )
            self.quitDialog()
            return

        overflow_right_flags, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.overflow_right_enabled_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for overflow right enabled!'
            )
            self.quitDialog()
            return

        poleni_right_values, ok = QgsVectorLayerUtils.getValues(input_layer, self.dialog.overflow_right_poleni_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for right poleni factor!'
            )
            self.quitDialog()
            return

        localbridgeheights, ok = QgsVectorLayerUtils.getValues(input_layer,
                                                               self.dialog.localbridgeheight_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for names!'
            )
            self.quitDialog()
            return

        bridgebodyheights, ok = QgsVectorLayerUtils.getValues(input_layer,
                                                              self.dialog.bridgebodyheight_box.expression(), os)
        if not ok:
            self.iface.messageBar().pushCritical(
                'River Profile Export',
                'Invalid expression for names!'
            )
            self.quitDialog()
            return

        dem_layer = self.dialog.raster_layer
        dem_band = self.dialog.raster_band_box.value()
        dem_name = self.dialog.raster_layer.name()
        dem_method = self.dialog.method_box.currentText()
        dem_nan = self.dialog.nan_box.value()
        dem_interpol = RasterInterpolator(dem_layer, dem_band, 1, 1, dem_method, dem_nan)
        dem_trans: object = QgsCoordinateTransform(input_layer.crs(), dem_layer.crs(), QgsProject.instance()).transform

        roughness_layer = self.dialog.roughness_layer
        roughness_band = self.dialog.roughness_band_box.value()
        roughness_nan = self.dialog.default_roughness_box.value()
        roughness_interpol = RasterInterpolator(roughness_layer, roughness_band, 10, 10, 'nearest neighbor (downscaling/upscaling)', roughness_nan)
        if roughness_layer:
            roughness_trans = QgsCoordinateTransform(input_layer.crs(), roughness_layer.crs(), QgsProject.instance())\
                .transform
        else:
            roughness_trans = lambda coord: coord

        progress = QProgressDialog('Exporting river profiles ...', 'Abort', 0, 2 * feature_count,
                                   self.iface.mainWindow())
        progress.setWindowTitle('River Profile Export')
        progress.canceled.connect(self.scheduleAbort)
        progress.show()

        text_blocks = {}
################################################################################
        #autostation
        if autostation:
            # linepre is the previous line
            linepre = []
            stationsum = 0
            autostations = list(range(0, len(profileids)))
            sortedindex=np.argsort(profileids)
            #bug fix
            feature = input_layer.getFeature(0)
            buff = feature.geometry()
            if not buff:
                layerinmemory=True
            else:
                layerinmemory = False
            for i, j in enumerate(sortedindex):
                if layerinmemory==True:
                    feature = input_layer.getFeature(j+1)
                else:
                    feature=input_layer.getFeature(j)
                line = []
                buff = feature.geometry()
                for p in buff.vertices():
                    line.append(p)
                #stationing starts with 0
                if i == 0:
                    autocalculatedstation = 0
                else:
                    #calculating the lowest hieght in the previous line
                    zpre = [dem_interpol(QgsPointXY(p)) for p in linepre]
                    #the location of the point with the lowest height
                    point_minpre = np.argmin(zpre)
                    znew = [dem_interpol(QgsPointXY(p)) for p in line]
                    point_minnew = np.argmin(znew)
                    difference = linepre[point_minpre].distance(line[point_minnew])
                    stationsum = difference + stationsum
                    autocalculatedstation = -stationsum
                linepre = line
                autostations[j]=autocalculatedstation

##################################################################################

        # iterate over profiles and extract attributes and points
        for i, f in enumerate(features):
            line = []
            buff = f.geometry()

            for p in buff.vertices():
                line.append(p)

            if flip_directions:
                line = list(reversed(line))

            if autostation:
                stations=autostations

            name = str(names[i])
            station = float(stations[i])
            # no check required
            delta_x = float(deltas[i])
            #check inflow, outflow, standard
            conn_type = str(conn_types[i]).lower()
            if check_river_prof_connection(conn_type) == 0:
                self.iface.messageBar().pushCritical(
                    'Error during river profile export',
                    'Invalid expression for river profile connection (standard/inflow/outflow) ({})! Aborting ...'
                    .format(conn_types[i])
                )
                self.quitDialog()
                return

            #check river, weir, bridge
            profile_type = str(profile_types[i]).lower()
            if check_river_prof_type(profile_type) == 0:
                self.iface.messageBar().pushCritical(
                    'Error during river profile export',
                    'Invalid expression for river profile type (river/bridge/weir) ({})! Aborting ...'
                    .format(profile_types[i])
                )
                self.quitDialog()
                return

            # no check required
            init_value = float(init_values[i])
            #check true, false
            point_bc = str(point_bc_flags[i]).lower()
            if check_true_false(point_bc) == 0:
                self.iface.messageBar().pushCritical(
                    'Error during river profile export',
                    'Invalid expression for point boundary condition flag (true/false) ({})! Aborting ...'
                    .format(point_bc_flags[i])
                )
                self.quitDialog()
                return
            # check true, false
            point_bc_stat = str(point_bc_stationary_flags[i]).lower()
            if check_true_false(point_bc_stat) ==0:
                self.iface.messageBar().pushCritical(
                    'Error during river profile export',
                    'Invalid expression for point boundary stationary flag (true/false) ({})! Aborting ...'
                    .format(point_bc_stationary_flags[i])
                )
                self.quitDialog()
                return
            # no check required
            point_bc_v = point_bc_values[i]
            # check true, false
            lat_bc = str(lateral_bc_flags[i]).lower()
            if check_true_false(lat_bc) == 0:
                self.iface.messageBar().pushCritical(
                    'Error during river profile export',
                    'Invalid expression for lateral boundary condition flag (true/false) ({})! Aborting ...'
                    .format(lateral_bc_flags[i])
                )
                self.quitDialog()
                return
            # check true, false
            lat_bc_stat = str(lateral_bc_stationary_flags[i]).lower()
            if check_true_false(lat_bc_stat) == 0:
                self.iface.messageBar().pushCritical(
                    'Error during river profile export',
                    'Invalid expression for lateral boundary stationary flag (true/false) ({})! Aborting ...'
                    .format(lateral_bc_stationary_flags[i])
                )
                self.quitDialog()
                return
            # no check required
            lat_bc_v = lateral_bc_values[i]
            # check true, false
            overflow_left = str(overflow_left_flags[i]).lower()
            if check_true_false(overflow_left) == 0:
                self.iface.messageBar().pushCritical(
                    'Error during river profile export',
                    'Invalid expression for left overflow flag (true/false) ({})! Aborting ...'
                    .format(overflow_left_flags[i])
                )
                self.quitDialog()
                return
            # no check required
            poleni_left = float(poleni_left_values[i])
            # check true, false
            overflow_right = str(overflow_right_flags[i]).lower()
            if check_true_false(overflow_right) == 0:
                self.iface.messageBar().pushCritical(
                    'Error during river profile export',
                    'Invalid expression for right overflow flag (true/false) ({})! Aborting ...'
                    .format(overflow_right_flags[i])
                )
                self.quitDialog()
                return
            # no check required
            poleni_right = float(poleni_right_values[i])

            # bridge value no check required; this are the individual values
            bridgebodyheight = float(bridgebodyheights[i])
            localbridgeheight = float(localbridgeheights[i])

            # label is None or empty
            if not name:
                self.iface.messageBar().pushCritical(
                    'Error during river profile export',
                    'Invalid name label found in field "{}"! Aborting ...'
                    .format(self.dialog.label_field_box.expression())
                )
                self.quitDialog()
                return

            # implicitly convert label to string
            name = str(name)
            # erase whitespace before
            name = name.replace(' ', '_')
            # label contains whitespaces
            if len(name.split(' ')) > 1:
                self.iface.messageBar().pushCritical(
                    'Error during river profile export',
                    'Labels must not contain whitespaces! Aborting ...'
                    .format(self.dialog.label_field_box.expression())
                )
                self.quitDialog()
                return

            if self.cancel:
                break

            # collect point data
            first = line[0]
            last = line[-1]

            x = [p.x() for p in line]
            y = [p.y() for p in line]
            z = [dem_interpol(QgsPointXY(p)) for p in line]
            mat = [int(roughness_interpol(roughness_trans(QgsPointXY(p)))) for p in line]

            if adjust_elevation:
                if z[0] <= z[1]:
                    z[0] = z[1] + 0.01
                if z[-1] <= z[-2]:
                    z[-1] = z[-2] + 0.01

            dist, ident = [], []
            # set channel type
            flag = False
            sumd=0

            #Main channel is used

            if self.dialog.channel_info == None:
                mainchannel_flag = False
            else:
                mainchannel_flag = True

                if self.dialog.channel_info.featureCount() != 1:
                    self.iface.messageBar().pushWarning('1-D River Profile Export',
                                                    'More than one polygon in main'
                                                    ' channel file; just the first one is used')

            for j, point in enumerate(line):
                if j==0:
                    d = 0
                else:
                    d = line[j-1].distance(point)
                sumd=d+sumd
                dist.append(sumd)

                if point == first:
                    ident.append(1)
                else:
                    # get main channel polygons
                    if mainchannel_flag == False:
                        ident.append(2)

                    else:
                        features_main = self.dialog.channel_info.getFeatures()
                        counter_p = 0
                        for poly in features_main:
                            counter_p = counter_p+1
                            if counter_p == 1:
                                geom_pol = poly.geometry()
                                if geom_pol.contains(QgsPointXY(point)) == True:
                                    ident.append(2)
                                    flag = True
                                elif flag == True:
                                    ident.append(3)
                                elif flag == False:
                                    ident.append(1)

                # elif point == last:
                #     ident.append(3)
                # else:
                #     ident.append(2)

            i_min = np.argmin(z)
            # minimum profile elevation
            z_min = z[int(i_min)]
            # init values are absolute values
            if abs_init:
                init_value -= z_min
            # write profile attributes
            block = ""
            block += 'ZONE T="{}" I={:d}\n'.format(name, len(line))
            block += 'AUXDATA ProfLDist="{:.2f}"\n'.format(station)
            block += 'AUXDATA DeltaXtable="{}"\n'.format(delta_x)
            if addfullriver:
                if stations[i]==max(stations):
                    block += 'AUXDATA ConnectionType="inflow"\n'
                elif stations[i]==min(stations):
                    block += 'AUXDATA ConnectionType="outflow"\n'
                else:
                    block += 'AUXDATA ConnectionType="{}"\n'.format(conn_type)
            else:
                block += 'AUXDATA ConnectionType="{}"\n'.format(conn_type)
            block += 'AUXDATA ProfType="{}"\n'.format(profile_type)
            block += 'AUXDATA InitCondition="{:.2f}"\n'.format(init_value)
            if addfullriver:
                if stations[i]==max(stations):
                    block += 'AUXDATA BoundaryPointCondition="true"\n'
                else:
                    block += 'AUXDATA BoundaryPointCondition="{}"\n'.format(str(point_bc).lower())
            else:
                block += 'AUXDATA BoundaryPointCondition="{}"\n'.format(str(point_bc).lower())
            block += 'AUXDATA BoundaryPointStationary="{}"\n'.format(str(point_bc_stat).lower())
            if addfullriver:
                if stations[i]==max(stations):
                    block += 'AUXDATA BoundaryPointValue="0.5"\n'
                else:
                    block += 'AUXDATA BoundaryPointValue="{}"\n'.format(point_bc_v)
            else:
                block += 'AUXDATA BoundaryPointValue="{}"\n'.format(point_bc_v)
            block += 'AUXDATA BoundaryLateralCondition="{}"\n'.format(str(lat_bc).lower())
            block += 'AUXDATA BoundaryLateralStationary="{}"\n'.format(str(lat_bc_stat).lower())
            block += 'AUXDATA BoundaryLateralValue="{}"\n'.format(lat_bc_v)
            if addfullriver or profile_types[i].lower()=="bridge":
                sortstations=np.sort(stations)
                if stations[i]==min(stations) or stations[i]==sortstations[1] or profile_types[i].lower()=="bridge":
                    block += 'AUXDATA OverflowCouplingLeft="false"\n'
                else:
                    block += 'AUXDATA OverflowCouplingLeft="{}"\n'.format(str(overflow_left).lower())
            else:
                block += 'AUXDATA OverflowCouplingLeft="{}"\n'.format(str(overflow_left).lower())
            block += 'AUXDATA PoleniFacLeft="{}"\n'.format(poleni_left)
            if addfullriver or profile_types[i].lower()=="bridge":
                sortstations = np.sort(stations)
                if stations[i]==min(stations) or stations[i]==sortstations[1] or profile_types[i].lower()=="bridge":
                    block += 'AUXDATA OverflowCouplingRight="false"\n'
                else:
                    block += 'AUXDATA OverflowCouplingRight="{}"\n'.format(str(overflow_right).lower())
            else:
                block += 'AUXDATA OverflowCouplingRight="{}"\n'.format(str(overflow_right).lower())
            if profile_types[i].lower()=="bridge":
                block += 'AUXDATA PoleniFacRight="{}"\n'.format(poleni_right)
                block += 'AUXDATA BridgeBodySize="{}"\n'.format(bridgebodyheight)
                block += 'AUXDATA LocalBridgeHeight="{}"\n\n'.format(localbridgeheight)
            else:
                block += 'AUXDATA PoleniFacRight="{}"\n\n'.format(poleni_right)

            # write profile points
            for xj, yj, zj, matj, distj, identj in zip(x, y, z, mat, dist, ident):
                block += '{x:.2f} {y:.2f} {z:.2f} {mat} {dist:.2f} {ident}\n'.format(x=xj, y=yj, z=zj, mat=matj,
                                                                                dist=distj, ident=identj)

            block += '\n'

            if station in text_blocks:
                self.iface.messageBar().pushWarning(
                    'River Profile Export',
                    'Duplicate profile station "{}"! Overwriting previous profile definition.'.format(station)
                )

            text_blocks[station] = block
            progress.setValue(i + 1)

            if self.cancel:
                break

        if self.cancel:
            progress.close()
            self.quitDialog()
            return

        # write (station-wise sorted) profiles and profile points to file
        with open(filename, 'w') as profile_file:
            profile_file.write('########################################################################\n')
            profile_file.write('# This file was automatically generated by ProMaiDes 1D-River Profile Export'
                                 '-QGIS-Plugin Version {version_1} \n'.format(version_1=VERSION))
            # date and time output
            now = datetime.now()
            dt_string = now.strftime("%d/%m/%Y %H:%M:%S")
            profile_file.write('# Generated at {dt_string_1} '.format(dt_string_1=dt_string))
            profile_file.write('from layer {filename_1} \n'.format(filename_1=input_layer.sourceName()))
            profile_file.write('#  based on height raster (DEM) {}  \n'.format(dem_name))
            if self.dialog.channel_info is not None:
                profile_file.write('#  based on main channel polygon {}  \n'.
                                   format(self.dialog.channel_info.sourceName()))
            if roughness_layer:
                profile_file.write('#  based on roughness raster {}  \n'.
                                   format(roughness_layer.name()))

            profile_file.write('# Comments are marked with #\n')
            profile_file.write('#\n')
            profile_file.write('# Explanation of data:\n')
            profile_file.write('#  The river profile-file is optimized for Tecplot; '
                                 ' each profile is one Zone with Name (T=..) NumberOfPoints (I=..) \n')
            profile_file.write('#  Several additional information to the profile follows; these '
                                 'attributes are valid to the distance to the next downstream profile \n')
            profile_file.write('#  The station value is from upstream to downstream decreasing\n')
            profile_file.write(
                '#  Finally the point data is given (always from left to right over profile in flow direction):\n '
                '#   x-coordinate y-coordinate z-coordinate material_id local_distance Lbank(1)_Main(2)_Rbank(3)\n'
                '#  Further remarks: \n#    Mterial_id and Lbank(1)_Main(2)_Rbank(3) are valid to '
                'the distance before the point\n#    First point has always Lbank(1) \n#    '
                'First point must be higher than point after; last point must be higher than point before\n#    '
                'If you want to have a higher point density, please use plugin Resample Polyline Vertices'
                ' \n')
            profile_file.write('#\n')
            profile_file.write('# Use in .ilm-file (just copy, delete the leading "#", set file(s)):\n')
            profile_file.write('#  In global section add 1 to !NOFRV = x+1 # Number of river models\n')
            profile_file.write('#  Set a river model between !$BEGINDESCRIPTION and !$ENDDESCRIPTION\n')
            profile_file.write('#   !$BEGINRVMODEL =  Index_(starts by 0) "NAME" \n')
            profile_file.write('#    !NOFPROF = 369 #Number of profiles in file (see below)\n')
            profile_file.write('#    !GEOMETRYFILE = "./PATH2FILE/FILE_NAME.txt"\n')
            profile_file.write('#    !INSTATBOUNDFILE = <SET>	\n')
            profile_file.write('#       $FILENAME="./PATH2FILE/FILE_NAME.txt" #Instationary boundary file of river\n')
            profile_file.write('#       $NUMBER_CURVE = 1 #number of curves in instationary boundary file\n')
            profile_file.write('#     </SET>	\n')
            profile_file.write('#    !LIMITS = <SET>	\n')
            profile_file.write('#       $RTOL = 1e-7 '
                               '#$RTOL = Defines relative tolerances [optional, standard value = 5.0e-7]\n')
            profile_file.write('#       $ATOL = 1e-6 #$ATOL = defines absolute tolerances '
                               '[optional, standard value = 5.0e-7]\n')
            profile_file.write('#$PROF_INTERFACE = false # false:= no interface is taken between '
                               'right/left bank and main channel; true:=interface is taken from the main channel '
                               '[optional, standard value = false]\n')
            profile_file.write('#$V_EXPLICIT = false # false:= velocity head is not taken into account for '
                               'calculation of the gradient; true:= it is taken into account'
                               ' [optional, standard value = false]	\n')
            profile_file.write('#     </SET>	\n')

            profile_file.write('#  !$ENDRVMODEL  \n')

            profile_file.write('########################################################################\n\n')

            profile_file.write('TITLE = "{}"\n'.format(input_layer.name()))
            profile_file.write('VARIABLES = "X", "Y", "Z", "MathType", "Distance", "Ident"\n')
            profile_file.write('DATASETAUXDATA NumOfProf = "{:d}"\n\n'.format(len(text_blocks)))

            # write profiles in reverse sorted order, i.e. decreasing station values
            for i, station in enumerate(reversed(sorted(text_blocks.keys()))):
                profile_file.write(text_blocks[station])
                progress.setValue(feature_count + i + 1)

        progress.close()

        #if not self.cancel:
        self.iface.messageBar().pushInfo('1D-River Profile Export', 'Export finished successfully!')

        self.quitDialog()


text:
import ROOT
import root_numpy
import numpy as np
from coffea.btag_tools import BTagScaleFactor

def ratio(h2_denom, h2_nom):

    denom = root_numpy.hist2array(h2_denom, return_edges=True)
    num = root_numpy.hist2array(h2_nom, return_edges=True)
    denom = np.nan_to_num(denom)
    num = np.nan_to_num(num)
    eff = np.divide(num[0], denom[0])
    eff = np.nan_to_num(eff)
    xedges, yedges = denom[1][0], denom[1][1]

    return eff, xedges, yedges

f_beff = ROOT.TFile("data/beff_test.root")
effs = {}
effs['b'] = f_beff.Get("b")
effs['c'] = f_beff.Get("c")
effs['usdg'] = f_beff.Get("usdg")

"""
def load_eff_maps(filename):

    f = ROOT.TFile(filename)
    h2_BTaggingEff_Denom_b = f.Get("h2_BTaggingEff_Denom_b")
    h2_BTaggingEff_Num_b = f.Get("h2_BTaggingEff_Num_b")
    h2_BTaggingEff_Denom_c = f.Get("h2_BTaggingEff_Denom_c")
    h2_BTaggingEff_Num_c = f.Get("h2_BTaggingEff_Num_c")
    h2_BTaggingEff_Denom_usdg = f.Get("h2_BTaggingEff_Denom_udsg")
    h2_BTaggingEff_Num_usdg = f.Get("h2_BTaggingEff_Num_udsg")

    effs = {}
    effs['b'], xedges, yedges = ratio(h2_BTaggingEff_Denom_b, h2_BTaggingEff_Num_b)
    effs['c'], xedges, yedges = ratio(h2_BTaggingEff_Denom_c, h2_BTaggingEff_Num_c)
    effs['usdg'], xedges, yedges = ratio(h2_BTaggingEff_Denom_usdg, h2_BTaggingEff_Num_usdg)

    return effs, xedges, yedges

effs, xedges, yedges = load_eff_maps('data/beff.root')
"""

btag_sf = BTagScaleFactor("data/CSV.csv", "medium", methods='comb,comb,comb')

def eval_sf_eff(ev):

    pt = ev["Jet_pt"]
    eta = ev["Jet_eta"]
    flavour = ev["Jet_flavour"]

    eff = []
    sf = []
    sf_up = []
    sf_down = []

    for i in range(len(pt)):

        if pt[i] > 1000:
            pt_clipped = 980.
        else:
            pt_clipped = pt[i]

        if abs(flavour[i])==5:
            flav='b'
            flav_id = 5
        elif abs(flavour[i])==4:
            flav='c'
            flav_id = 4
        else:
            flav='usdg'
            flav_id = 0

        """
        ind_pt = np.digitize(pt_clipped, xedges)
        ind_eta = np.digitize(eta[i], yedges)
        """

        try:
            #eff.append( effs[flav][ind_pt, ind_eta] )
            eff.append( effs[flav].GetBinContent( effs[flav].GetXaxis().FindBin(pt_clipped),
                                                  effs[flav].GetYaxis().FindBin(eta[i])) )
            sf.append(btag_sf.eval("central", flav_id, abs(eta[i]), pt_clipped, ignore_missing=True))
            sf_up.append(btag_sf.eval("up", flav_id, abs(eta[i]), pt_clipped, ignore_missing=True))
            sf_down.append(btag_sf.eval("down", flav_id, abs(eta[i]), pt_clipped, ignore_missing=True))
        except:
            print( "Error in calculating the b-tagging SF:", pt_clipped, eta[i], flavour[i] )

    return {"Jet_btagSF" : sf, "Jet_btagSF_up" : sf_up, "Jet_btagSF_down" : sf_down, "Jet_beff" : eff}

"""
def eval_btagSF(ev):

    # !! CHECK FLAVOUR
    pt = ev["Jet_pt"]
    eta = ev["Jet_eta"]
    flavour = ev["Jet_flavour"]

    sf = btag_sf.eval("central", flavour, eta, pt, ignore_missing=True)

    return sf
"""

# Corresponds to https://twiki.cern.ch/twiki/bin/viewauth/CMS/BTagSFMethods 1A
def b_weight_method1(ev, syst='cent', njets = -1):

    jet_eff = ev["Jet_beff"]
    if syst == 'cent':
        jet_btagSF = ev["Jet_btagSF"]
    elif syst == 'up':
        jet_btagSF = ev["Jet_btagSF_up"]
    elif syst == 'down':
        jet_btagSF = ev["Jet_btagSF_down"]
    else:
        print("No valid b-tag variation!!")
    jet_csvDisc = ev["Jet_csvDisc"]

    pMC=1.0
    pData=1.0
    weight=1.0

    iJet = 0
    for eff, csvDisc, sf in zip(jet_eff, jet_csvDisc, jet_btagSF):

        if iJet == njets:
            break

        if csvDisc > 0.679:
            pMC=pMC*eff
            pData=pData*(sf*eff)

        else:
            pMC=pMC*(1.0 - eff)
            pData=pData*(1.0 - sf*eff)

        iJet += 1

    if pMC != 0.0:
        weight = pData/pMC

    return weight

# Corresponds to https://twiki.cern.ch/twiki/bin/viewauth/CMS/BTagSFMethods 1B
def b_weight_method2(ev, syst='cent', njets = -1):

    jet_eff = ev["Jet_beff"]
    if syst == 'cent':
        jet_btagSF = ev["Jet_btagSF"]
    elif syst == 'up':
        jet_btagSF = ev["Jet_btagSF_up"]
    elif syst == 'down':
        jet_btagSF = ev["Jet_btagSF_down"]
    else:
        print("No valid b-tag variation!!")
    jet_csvDisc = ev["Jet_csvDisc"]

    weight_0tags = 1.0

    iJet = 0
    for eff, csvDisc, sf in zip(jet_eff, jet_csvDisc, jet_btagSF):

        if iJet == njets:
            break
        try:
            weight_0tags = weight_0tags * (1 - sf * eff)
        except:
            print("sf", sf, "eff", eff)

        iJet += 1

    weight_at_least_1tag = 1.0 - weight_0tags

    return weight_at_least_1tag

def count_btags(jet_csvDisc, njets = -1):

    if njets < 0:
        return sum(jet_csvDisc>0.679)
    else:
        return sum(jet_csvDisc[0:njets] > 0.679)

def at_least_1tag(df):
    return df[df["Jet_nbtags"]>0]

def no_tag(df):
    return df[df["Jet_nbtags"]<1]

ext: py

max_stars_repo_name: llayer/cmsopen

max_issues_repo_name: llayer/cmsopen

max_forks_repo_name: llayer/cmsopen

reason: import numpy


text:
import argparse
import csv
import os
import numpy as np
import scipy as sc
import loompy

def generate_col_attr(qc_paths):
    """Converts the QC of Smart Seq2 gene file pipeline outputs to loom file
    Args:
        qc_path (str): path to the QCs csv
    """
    # read the QC values
    qc_path = [p for p in qc_paths if p.endswith("_QCs.csv")][0]
    with open(qc_path, 'r') as f:
        qc_values = [row for row in csv.reader(f)]

    metadata_labels = qc_values[0][1:]
    metadata_values = qc_values[2][1:]
    cell_id = qc_values[2][0]

    string_metadata = {}
    numeric_metadata = {}

    for label, value in zip(metadata_labels, metadata_values):

        # See if this is numeric or string
        numeric_value = None
        try:
            numeric_value = float(value)
        except ValueError:
            try:
                numeric_value = float(value.strip("%"))/100
            except ValueError:
                pass

        if numeric_value is not None:
            numeric_metadata[label] = numeric_value
        else:
            string_metadata[label] = value

    # Metrics
    # Write the string and numeric metadata separately
    sorted_string_labels = sorted(string_metadata.keys())
    sorted_string_values = [string_metadata[m] for m in sorted_string_labels]
    sorted_numeric_labels = sorted(numeric_metadata.keys())
    sorted_numeric_values = [numeric_metadata[m] for m in sorted_numeric_labels]

    # Column attributes
    col_attrs = dict()
    col_attrs["cell_names"] = [cell_id]

    numeric_field_names = np.array(sorted_numeric_labels[:])
    for i in range(0, numeric_field_names.shape[0]):
        name = numeric_field_names[i]
        data = np.array([sorted_numeric_values])[:,i]
        col_attrs[name] = data
    string_field_names = np.array(sorted_string_labels)
    for i in range(0, string_field_names.shape[0]):
        name = string_field_names[i]
        data = np.array([sorted_string_values])[:,i]
        col_attrs[name] = data

    return col_attrs

def generate_csr_spase_coo(expr):

    nrows, ncols = np.shape(expr)
    expr_coo = sc.sparse.coo_matrix(expr[:])
    xcoord = []
    ycoord = []
    value = []

    for k in range(0, expr_coo.data.shape[0]):
        xcoord.append(expr_coo.row[k])
        ycoord.append(expr_coo.col[k])
        value.append(expr_coo.data[k])

    xcoord = np.asarray(xcoord)
    ycoord = np.asarray(ycoord)
    value = np.asarray(value)

    expr_sp_t = sc.sparse.coo_matrix((value, (ycoord, xcoord)), shape=(ncols, nrows))
    return expr_sp_t

def generate_row_attr_and_matrix(rsem_gene_results_path):
    """Converts the Smart Seq2 gene file pipeline outputs to loom file
    Args:
        input_path (str): file where the SS2 pipeline expression counts are
    """
    reader = csv.DictReader(open(rsem_gene_results_path), delimiter="\t")
    expression_values = {}
    count_values = {}
    for row in reader:
        expression_values[row["gene_id"]] = float(row["TPM"])
        count_values[row["gene_id"]] = float(row["expected_count"])
    sorted_gene_ids = sorted(expression_values.keys())
    sorted_tpms = [expression_values[g] for g in sorted_gene_ids]
    sorted_counts = [count_values[g] for g in sorted_gene_ids]
    expression_tpms = generate_csr_spase_coo( [sorted_tpms])
    expected_counts = generate_csr_spase_coo([sorted_counts])
    row_attrs = {"ensembl_ids":np.array(sorted_gene_ids),"gene_names":np.array(sorted_gene_ids)}

    return row_attrs, expression_tpms,expected_counts

def create_loom_files(sample_id, qc_files, rsem_genes_results_file,
                      output_loom_path):
    """This function creates the loom file or folder structure in output_loom_path in
       format file_format, with sample_id from the input folder analysis_output_path
    Args:
        sample_id (str): sample or cell id
        qc_analysis_output_files_string (str): a string with the file names in the QCGroup of SS2
            pipeline output, separated by commas
        rsem_genes_results_file (str): the file for the expression count
        output_loom_path (str): location of the output loom
    """
    # generate a dictionarty of column attributes
    col_attrs =  generate_col_attr(qc_files)

    # add the expression count matrix data
    # generate a dictionary of row attributes
    row_attrs, expr_tpms, expr_counts = generate_row_attr_and_matrix(rsem_genes_results_file)

    attrDict = dict()
    attrDict['sample_id'] = sample_id

    #generate loom file
    loompy.create(output_loom_path, expr_tpms, row_attrs, col_attrs, file_attrs=attrDict)
    ds = loompy.connect(output_loom_path)
    ds.layers['estimated_counts'] = expr_counts
    ds.close()

def main():
    description = """This script converts the some of the SmartSeq2 pipeline outputs in to
                   loom format (https://linnarssonlab.org/loompy/format/index.html) relevant output.
                   This script can be used as a module or run as a command line script."""

    parser = argparse.ArgumentParser(description=description)
    parser.add_argument('--qc_files',
                        nargs = "+",
                        help=('the grouped QC files from the GroupQCOutputs task of SS2 '
                              'Single Sample workflow'))

    parser.add_argument('--rsem_genes_results',
                        help='path to the folder containing the files to be added to the loom')

    parser.add_argument('--output_loom_path',
                        help='path where the loom file is to be created')

    parser.add_argument('--sample_id',
                        default="Unknown sample",
                        help='the sample name in the bundle')

    args = parser.parse_args()

    create_loom_files(args.sample_id, args.qc_files, args.rsem_genes_results, args.output_loom_path)

if __name__ == '__main__':
    main()

max_stars_repo_name: truwl/skylab

max_issues_repo_name: truwl/skylab

max_forks_repo_name: truwl/skylab

reason: import numpy,import scipy


text:
#!/usr/bin/env python

"""Basic usage of the pybotics package."""
import numpy as np
from mdh import __version__ as version
from mdh.link import mdh_params, JointType
from mdh.kinematic_chain import KinematicChain

pi = np.pi

print(f">> {version}")

# make it print better
np.set_printoptions(suppress=True)

def kevin():
    mm = pi
    """alpha a theta d"""
    dh = [
        {'alpha':    0, 'a':  0, 'theta': 0, 'd': 0, 'type': 1, "max_min": [mm,-mm]},
        {'alpha': pi/2, 'a': 52, 'theta': 0, 'd': 0, 'type': 1, "max_min": [mm,-mm]},
        {'alpha':    0, 'a': 89, 'theta': 0, 'd': 0, 'type': 1, "max_min": [mm,-mm]},
        {'alpha':    0, 'a': 90, 'theta': 0, 'd': 0, 'type': 1, "max_min": [mm,-mm]},
        {'alpha':    0, 'a': 95, 'theta': 0, 'd': 0, 'type': 1, "max_min": [mm,-mm]}
    ]

    return dh

def main():

    t = kevin()
    kc = KinematicChain.from_parameters(t)

    # j = [
    #     [0,0,0,0,0],
    #     [0,pi/4,-pi/4,-pi/2,0], # 204, 0, -32
    #     # [0,0,-pi/2,0,0],
    #     # np.deg2rad([90,131.17,-106.19,-114.98]),
    #     np.deg2rad([-45.00, 77.41, -98.15, -69.27, 0]) # 110, -110, -40
    # ]

    # for jj in j:
    #     # print(jj)
    #     t = kc.transform(jj)
    #     print(t)
    #     print("----------------------")

    j = [
        [204.93, 0, -32.067],
        [110,-110,-40],
        [110,110,-40],
        [110,0,-70]
    ]

    for jj in j:
        p = np.array([
            [0,1,0, jj[0]],
            [0,0,-1,jj[1]],
            [-1,0,0,jj[2]],
            [0,0,0,1]
        ])
        rads = kc.inverse(p)
        print(">>", np.rad2deg(rads))

if __name__ == "__main__":
    main()

max_stars_repo_name: MultipedRobotics/dh

max_issues_repo_name: MultipedRobotics/dh

max_forks_repo_name: MultipedRobotics/dh


text:
import cv2
import numpy as np

def abs_sobel_thresh(img, orient='x', thresh=(0, 255),sobel_kernel=3,):
#     gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    sobel = None
    if(orient == 'x'):
        sobel = cv2.Sobel(img, cv2.CV_64F, 1, 0,ksize=sobel_kernel)
    else:
        sobel = cv2.Sobel(img, cv2.CV_64F, 0, 1,ksize=sobel_kernel)
    abs_sobel = np.absolute(sobel)
    scaled_sobel = np.uint8(255*abs_sobel/np.max(abs_sobel))
    sxbinary = np.zeros_like(scaled_sobel)
    sxbinary[(scaled_sobel >= thresh[0]) & (scaled_sobel <= thresh[1])] = 1
    return sxbinary

def mag_thresh(img, sobel_kernel=3, mag_thresh=(0, 255)):

#     gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    sobelx = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=sobel_kernel)
    sobely = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=sobel_kernel)
    magnitude = np.sqrt(sobelx**2 + sobely**2)

    scaled_sobel = np.uint8(255*magnitude/np.max(magnitude))

    binary = np.zeros_like(scaled_sobel)
    binary[(scaled_sobel >= mag_thresh[0]) & (scaled_sobel <= mag_thresh[1])] = 1
    return binary

def dir_threshold(img, sobel_kernel=3, thresh=(0, np.pi/2)):
#     gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    sobelx = np.absolute(cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=sobel_kernel))
    sobely = np.absolute(cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=sobel_kernel))
    direction = np.arctan2(sobely, sobelx)

    binary = np.zeros_like(direction)
    binary[(direction >= thresh[0]) & (direction <= thresh[1])] = 1
    return binary

def filter_with_hough(binary):
    output = np.zeros_like(binary)
    hough = cv2.HoughLinesP(binary, 1, (np.pi/180)*1, 1, np.array([]), minLineLength=8, maxLineGap=5)
    for x1,y1,x2,y2 in [x[0] for x in hough]:
        cv2.line(output,(x1,y1),(x2,y2),(255,255,255),2)
    return output

def color_treshold(binary, treshold=(170,255)):
    thresh_min = treshold[0]
    thresh_max = treshold[1]
    r_binary = np.zeros_like(binary)
    r_binary[(binary >= thresh_min) & (binary <= thresh_max)] = 1
    return r_binary

def get_tresholded_img(image):
    ksize = 15
    gkernel_size = 17
    image = cv2.GaussianBlur(image, (gkernel_size, gkernel_size), 0)
#     return image
    hls = cv2.cvtColor(image, cv2.COLOR_RGB2HLS)
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    lab = cv2.cvtColor(image, cv2.COLOR_BGR2Lab)
    luv = cv2.cvtColor(image, cv2.COLOR_BGR2LUV)

#     L = lab[:,:,0]
#     a = lab[:,:,1]
    lab_b = lab[:,:,2]
    luv_l = luv[:,:,0]

    l_channel = hls[:,:,1]
    s_channel = hls[:,:,2]
#     h_channel = hls[:,:,0]
    l_channel = cv2.equalizeHist(l_channel)
    s_channel = cv2.equalizeHist(s_channel)
    lab_b = cv2.equalizeHist(lab_b)
    luv_l = cv2.equalizeHist(luv_l)
#     h_channel = cv2.equalizeHist(s_channel)
    gray = cv2.equalizeHist(gray)

#     L = cv2.equalizeHist(L)

    working_channel = l_channel + s_channel
#     mag_binary = mag_thresh(s_channel, sobel_kernel=ksize, mag_thresh=(25, 100))
#     s_binary = color_treshold(s_channel, (170,255))
#     gradx = abs_sobel_thresh(working_channel, orient='x', sobel_kernel=ksize, thresh=(25, 100))
#     gradx = abs_sobel_thresh(s_channel, orient='x', sobel_kernel=ksize, thresh=(25, 100))
#     grady = abs_sobel_thresh(working_channel, orient='y', sobel_kernel=ksize, thresh=(20, 100))
    dir_binary = dir_threshold(s_channel, sobel_kernel=ksize, thresh=(0.4, 1.3))
    gradx = abs_sobel_thresh(working_channel, orient='x', sobel_kernel=ksize, thresh=(25, 100))

    lab_b = color_treshold(lab_b, (249,255))
    luv_l = color_treshold(luv_l, (249,255))
    hls_s = color_treshold(s_channel, (249,255))
    combined = np.zeros_like(s_channel)
#     | (gradx*dir_binary == 1)
    combined[(lab_b == 1)|(luv_l ==1)| (hls_s == 1)] = 1
    combined_image = np.array(combined * 255, dtype = np.uint8)
    return combined_image


text:
import os
import sys

import cv2
import numpy as np

class TumDataset:
    def __init__(self, path):
        rgb_path, depth_path = self.provide_rgb_and_depth_path(path)

        rgb_filenames, depth_filenames = self.provide_filenames(rgb_path, depth_path)
        rgb_filenames.sort()
        depth_filenames.sort()

        self.depth_images = [os.path.join(depth_path, filename) for filename in depth_filenames]
        self.rgb_images = [os.path.join(rgb_path, filename) for filename in rgb_filenames]
        self.depth_to_rgb_index = []

        self.match_rgb_with_depth(rgb_filenames, depth_filenames)

    def provide_rgb_and_depth_path(self, path):
        depth_path = os.path.join(path, "depth")
        rgb_path = os.path.join(path, "rgb")

        return rgb_path, depth_path

    def provide_filenames(self, rgb_path, depth_path):
        rgb_filenames = os.listdir(rgb_path)
        depth_filenames = os.listdir(depth_path)

        return rgb_filenames, depth_filenames

    def read_depth_image(self, frame_num):
        depth_frame_path = self.depth_images[frame_num]
        return cv2.imread(depth_frame_path, cv2.IMREAD_ANYDEPTH)

    def match_rgb_with_depth(self, rgb_filenames, depth_filenames):
        rgb_index = 0
        depth_index = 0
        prev_delta = float('inf')
        while depth_index < len(depth_filenames) and rgb_index < len(rgb_filenames):
            rgb_timestamp = float(rgb_filenames[rgb_index][:-4])
            depth_timestamp = float(depth_filenames[depth_index][:-4])
            delta = abs(depth_timestamp - rgb_timestamp)

            if rgb_timestamp <= depth_timestamp:
                prev_delta = delta
                rgb_index += 1
                continue

            if prev_delta < delta:
                self.depth_to_rgb_index.append(rgb_index - 1)
            else:
                self.depth_to_rgb_index.append(rgb_index)

            depth_index += 1

class CameraIntrinsics:
    def __init__(self, width, height, fx, fy, cx, cy):
        self.width = width
        self.height = height
        self.fx = fx
        self.fy = fy
        self.cx = cx
        self.cy = cy

def get_cam_intrinsics(image_shape):
    return CameraIntrinsics(
            width=image_shape[1],
            height=image_shape[0],
            fx=591.1,  # X-axis focal length
            fy=590.1,  # Y-axis focal length
            cx=331.0,  # X-axis principle point
            cy=234.0  # Y-axis principle point
        )

if __name__ == "__main__":
    path = sys.argv[1]
    tum_loader = TumDataset(path)
    path_out = os.path.join(path, "out")
    if not os.path.exists(path_out):
        os.mkdir(path_out)

    i = 0
    for index, depth_path in enumerate(tum_loader.depth_images):
        rgb_index = tum_loader.depth_to_rgb_index[index]
        rgb_path = tum_loader.rgb_images[rgb_index]

        rgb_image = cv2.imread(rgb_path)
        depth_image = tum_loader.read_depth_image(index)
        if index == 0:
            camera_intrinsics = get_cam_intrinsics(depth_image.shape)
            column_indices = np.tile(np.arange(camera_intrinsics.width), (camera_intrinsics.height, 1)).flatten()
            row_indices = np.transpose(np.tile(np.arange(camera_intrinsics.height), (camera_intrinsics.width, 1))).flatten()

        rgb_image[depth_image == 0] = [0, 0, 0]

        factor = 5000
        points = np.zeros((camera_intrinsics.width * camera_intrinsics.height, 3))

        points[:, 2] = depth_image.flatten() / factor
        points[:, 0] = (column_indices - camera_intrinsics.cx) * points[:, 2] / camera_intrinsics.fx
        points[:, 1] = (row_indices - camera_intrinsics.cy) * points[:, 2] / camera_intrinsics.fy
        distance = np.sqrt(np.sum(np.square(points), axis=1)).reshape((480, 640))
        rgb_image[distance >= 5] = [0, 0, 0]

        cv2.imwrite(os.path.join(path_out, "{:05d}.png".format(i)), rgb_image)
        i += 1


text:
import nn
import numpy as np
import sys

from util import *
from visualize import *
from layers import *

# XTrain - List of training input Data
# YTrain - Corresponding list of training data labels
# XVal - List of validation input Data
# YVal - Corresponding list of validation data labels
# XTest - List of testing input Data
# YTest - Corresponding list of testing data labels

def taskSquare(draw):
	XTrain, YTrain, XVal, YVal, XTest, YTest = readSquare()
	# Create a NeuralNetwork object 'nn1' as follows with optimal parameters. For parameter definition, refer to nn.py file.
	# nn1 = nn.NeuralNetwork(out_nodes, alpha, batchSize, epochs)
	# Add layers to neural network corresponding to inputs and outputs of given data
	# Eg. nn1.addLayer(FullyConnectedLayer(x,y))
	###############################################
	# TASK 2.1 - YOUR CODE HERE
	nn1 = nn.NeuralNetwork(2, 0.3, 4, 20)
	nn1.addLayer(FullyConnectedLayer(2,4))
	nn1.addLayer(FullyConnectedLayer(4,2))
	###############################################
	nn1.train(XTrain, YTrain, XVal, YVal, False, True)
	pred, acc = nn1.validate(XTest, YTest)
	print('Test Accuracy ',acc)
	# Run script visualizeTruth.py to visualize ground truth. Run command 'python3 visualizeTruth.py 2'
	# Use drawSquare(XTest, pred) to visualize YOUR predictions.
	# if draw:
	# 	drawSquare(XTest, pred)
	return nn1, XTest, YTest

def taskSemiCircle(draw):
	XTrain, YTrain, XVal, YVal, XTest, YTest = readSemiCircle()
	# Create a NeuralNetwork object 'nn1' as follows with optimal parameters. For parameter definition, refer to nn.py file.
	# nn1 = nn.NeuralNetwork(out_nodes, alpha, batchSize, epochs)
	# Add layers to neural network corresponding to inputs and outputs of given data
	# Eg. nn1.addLayer(FullyConnectedLayer(x,y))
	###############################################
	# TASK 2.2 - YOUR CODE HERE
	nn1 = nn.NeuralNetwork(2, 0.3, 4, 20)
	nn1.addLayer(FullyConnectedLayer(2,2))
	nn1.addLayer(FullyConnectedLayer(2,2))
	###############################################
	nn1.train(XTrain, YTrain, XVal, YVal, False, True)
	pred, acc  = nn1.validate(XTest, YTest)
	print('Test Accuracy ',acc)
	# Run script visualizeTruth.py to visualize ground truth. Run command 'python3 visualizeTruth.py 4'
	# Use drawSemiCircle(XTest, pred) to vnisualize YOUR predictions.
	if draw:
		drawSemiCircle(XTest, pred)
	return nn1, XTest, YTest

def taskMnist():
	XTrain, YTrain, XVal, YVal, XTest, YTest = readMNIST()
	# Create a NeuralNetwork object 'nn1' as follows with optimal parameters. For parameter definition, refer to nn.py file.
	# nn1 = nn.NeuralNetwork(out_nodes, alpha, batchSize, epochs)
	# Add layers to neural network corresponding to inputs and outputs of given data
	# Eg. nn1.addLayer(FullyConnectedLayer(x,y))
	###############################################
	# TASK 2.3 - YOUR CODE HERE
	# seed value = 30
	nn1 = nn.NeuralNetwork(10, 0.3, 4, 20)
	nn1.addLayer(FullyConnectedLayer(784,40))
	nn1.addLayer(FullyConnectedLayer(40,10))
	###############################################
	nn1.train(XTrain, YTrain, XVal, YVal, False, True)
	pred, acc  = nn1.validate(XTest, YTest)
	print('Test Accuracy ',acc)
	return nn1, XTest, YTest

def taskCifar10():
	XTrain, YTrain, XVal, YVal, XTest, YTest = readCIFAR10()

	XTrain = XTrain[0:5000,:,:,:]
	XVal = XVal[0:1000,:,:,:]
	XTest = XTest[0:1000,:,:,:]
	YVal = YVal[0:1000,:]
	YTest = YTest[0:1000,:]
	YTrain = YTrain[0:5000,:]

	modelName = 'model.npy'
	# # Create a NeuralNetwork object 'nn1' as follows with optimal parameters. For parameter definition, refer to nn.py file.
	# # nn1 = nn.NeuralNetwork(out_nodes, alpha, batchSize, epochs)
	# # Add layers to neural network corresponding to inputs and outputs of given data
	# # Eg. nn1.addLayer(FullyConnectedLayer(x,y))
	# ###############################################
	# # TASK 2.4 - YOUR CODE HERE
	nn1 = nn.NeuralNetwork(10,0.25,32,25)
	nn1.addLayer(ConvolutionLayer([3,32,32],[10,10],10,3))
	nn1.addLayer(AvgPoolingLayer([10,8,8],[2,2], 2))
	nn1.addLayer(ConvolutionLayer([10,4,4],[2,2],15,2))
	nn1.addLayer(FlattenLayer())
	nn1.addLayer(FullyConnectedLayer(60,10))
	###################################################
	return nn1,  XTest, YTest, modelName # UNCOMMENT THIS LINE WHILE SUBMISSION

	nn1.train(XTrain, YTrain, XVal, YVal, True, True, loadModel=False, saveModel=True, modelName=modelName)
	pred, acc = nn1.validate(XTest, YTest)
	print('Test Accuracy ',acc)


text:
#!/usr/bin/python
# coding: utf8
'''
Created on Oct 27, 2010
Update  on 2017-05-18
Logistic Regression Working Module
Author: Peter Harrington/羊三/小瑶
GitHub: https://github.com/apachecn/MachineLearning
'''
from numpy import *
import matplotlib.pyplot as plt

# ---------------------------------------------------------------------------
# 使用 Logistic 回归在简单数据集上的分类

# 解析数据
def loadDataSet(file_name):
    '''
    Desc:
        加载并解析数据
    Args:
        file_name -- 文件名称，要解析的文件所在磁盘位置
    Returns:
        dataMat -- 原始数据的特征
        labelMat -- 原始数据的标签，也就是每条样本对应的类别
    '''
    # dataMat为原始数据， labelMat为原始数据的标签
    dataMat = []
    labelMat = []
    fr = open(file_name)
    for line in fr.readlines():
        lineArr = line.strip().split()
        if len(lineArr) == 1:
            continue    # 这里如果就一个空的元素，则跳过本次循环
        # 为了方便计算，我们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 作为 X0
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

# sigmoid跳跃函数
def sigmoid(inX):
    # return 1.0 / (1 + exp(-inX))

    # Tanh是Sigmoid的变形，与 sigmoid 不同的是，tanh 是0均值的。因此，实际应用中，tanh 会比 sigmoid 更好。
    return 2 * 1.0/(1+exp(-2*inX)) - 1

# 正常的处理方案
# 两个参数：第一个参数==> dataMatIn 是一个2维NumPy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。
# 第二个参数==> classLabels 是类别标签，它是一个 1*100 的行向量。为了便于矩阵计算，需要将该行向量转换为列向量，做法是将原向量转置，再将它赋值给labelMat。
def gradAscent(dataMatIn, classLabels):
    '''
    Desc:
        正常的梯度上升法
    Args:
        dataMatIn -- 输入的 数据的特征 List
        classLabels -- 输入的数据的类别标签
    Returns:
        array(weights) -- 得到的最佳回归系数
    '''

    # 转化为矩阵[[1,1,2],[1,1,2]....]
    dataMatrix = mat(dataMatIn)  # 转换为 NumPy 矩阵
    # 转化为矩阵[[0,1,0,1,0,1.....]]，并转制[[0],[1],[0].....]
    # transpose() 行列转置函数
    # 将行向量转化为列向量   =>  矩阵的转置
    labelMat = mat(classLabels).transpose()  # 首先将数组转换为 NumPy 矩阵，然后再将行向量转置为列向量
    # m->数据量，样本数 n->特征数
    m, n = shape(dataMatrix)
    # print m, n, '__'*10, shape(dataMatrix.transpose()), '__'*100
    # alpha代表向目标移动的步长
    alpha = 0.001
    # 迭代次数
    maxCycles = 500
    # 生成一个长度和特征数相同的矩阵，此处n为3 -> [[1],[1],[1]]
    # weights 代表回归系数， 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵，其中的数全部都是 1
    weights = ones((n, 1))
    for k in range(maxCycles):  # heavy on matrix operations
        # m*3 的矩阵 * 3*1 的单位矩阵 ＝ m*1的矩阵
        # 那么乘上单位矩阵的意义，就代表：通过公式得到的理论值
        # 参考地址： 矩阵乘法的本质是什么？ https://www.zhihu.com/question/21351965/answer/31050145
        # print 'dataMatrix====', dataMatrix
        # print 'weights====', weights
        # n*3   *  3*1  = n*1
        h = sigmoid(dataMatrix * weights)  # 矩阵乘法
        # print 'hhhhhhh====', h
        # labelMat是实际值
        error = (labelMat - h)  # 向量相减
        # 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况，最后得出 x1,x2,xn的系数的偏移量
        weights = weights + alpha * dataMatrix.transpose() * error  # 矩阵乘法，最后得到回归系数
    return array(weights)

# 随机梯度下降
# 梯度下降优化算法在每次更新数据集时都需要遍历整个数据集，计算复杂都较高
# 随机梯度下降一次只用一个样本点来更新回归系数
def stocGradAscent0(dataMatrix, classLabels):
    '''
    Desc:
        随机梯度下降，只使用一个样本点来更新回归系数
    Args:
        dataMatrix -- 输入数据的数据特征（除去最后一列）
        classLabels -- 输入数据的类别标签（最后一列数据）
    Returns:
        weights -- 得到的最佳回归系数
    '''
    m, n = shape(dataMatrix)
    alpha = 0.01
    # n*1的矩阵
    # 函数ones创建一个全1的数组
    weights = ones(n)  # 初始化长度为n的数组，元素全部为 1
    for i in range(m):
        # sum(dataMatrix[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn,此处求出的 h 是一个具体的数值，而不是一个矩阵
        h = sigmoid(sum(dataMatrix[i] * weights))
        # print 'dataMatrix[i]===', dataMatrix[i]
        # 计算真实类别与预测类别之间的差值，然后按照该差值调整回归系数
        error = classLabels[i] - h
        # 0.01*(1*1)*(1*n)
        # print weights, "*" * 10, dataMatrix[i], "*" * 10, error
        weights = weights + alpha * error * dataMatrix[i]
    return weights

# 随机梯度下降算法（随机化）
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    '''
    Desc:
        改进版的随机梯度下降，使用随机的一个样本来更新回归系数
    Args:
        dataMatrix -- 输入数据的数据特征（除去最后一列数据）
        classLabels -- 输入数据的类别标签（最后一列数据）
        numIter=150 --  迭代次数
    Returns:
        weights -- 得到的最佳回归系数
    '''
    m, n = shape(dataMatrix)
    weights = ones(n)  # 创建与列数相同的矩阵的系数矩阵，所有的元素都是1
    # 随机梯度, 循环150,观察是否收敛
    for j in range(numIter):
        # [0, 1, 2 .. m-1]
        dataIndex = range(m)
        for i in range(m):
            # i和j的不断增大，导致alpha的值不断减少，但是不为0
            alpha = 4 / (
                1.0 + j + i
            ) + 0.0001  # alpha 会随着迭代不断减小，但永远不会减小到0，因为后边还有一个常数项0.0001
            # 随机产生一个 0～len()之间的一个值
            # random.uniform(x, y) 方法将随机生成下一个实数，它在[x,y]范围内,x是这个范围内的最小值，y是这个范围内的最大值。
            randIndex = int(random.uniform(0, len(dataIndex)))
            # sum(dataMatrix[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn
            h = sigmoid(sum(dataMatrix[dataIndex[randIndex]] * weights))
            error = classLabels[dataIndex[randIndex]] - h
            # print weights, '__h=%s' % h, '__'*20, alpha, '__'*20, error, '__'*20, dataMatrix[randIndex]
            weights = weights + alpha * error * dataMatrix[dataIndex[randIndex]]
            del (dataIndex[randIndex])
    return weights

# 可视化展示
def plotBestFit(dataArr, labelMat, weights):
    '''
        Desc:
            将我们得到的数据可视化展示出来
        Args:
            dataArr:样本数据的特征
            labelMat:样本数据的类别标签，即目标变量
            weights:回归系数
        Returns:
            None
    '''

    n = shape(dataArr)[0]
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):
        if int(labelMat[i]) == 1:
            xcord1.append(dataArr[i, 1])
            ycord1.append(dataArr[i, 2])
        else:
            xcord2.append(dataArr[i, 1])
            ycord2.append(dataArr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    x = arange(-3.0, 3.0, 0.1)
    """
    y的由来，卧槽，是不是没看懂？
    首先理论上是这个样子的。
    dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
    w0*x0+w1*x1+w2*x2=f(x)
    x0最开始就设置为1叻， x2就是我们画图的y值，而f(x)被我们磨合误差给算到w0,w1,w2身上去了
    所以： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2
    """
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y)
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.show()

def simpleTest():
    # 1.收集并准备数据
    dataMat, labelMat = loadDataSet("input/5.Logistic/TestSet.txt")

    # print dataMat, '---\n', labelMat
    # 2.训练模型，  f(x)=a1*x1+b2*x2+..+nn*xn中 (a1,b2, .., nn).T的矩阵值
    # 因为数组没有是复制n份， array的乘法就是乘法
    dataArr = array(dataMat)
    # print dataArr
    # weights = gradAscent(dataArr, labelMat)
    # weights = stocGradAscent0(dataArr, labelMat)
    weights = stocGradAscent1(dataArr, labelMat)
    # print '*'*30, weights

    # 数据可视化
    plotBestFit(dataArr, labelMat, weights)

# --------------------------------------------------------------------------------
# 从疝气病症预测病马的死亡率
# 分类函数，根据回归系数和特征向量来计算 Sigmoid的值
def classifyVector(inX, weights):
    '''
    Desc:
        最终的分类函数，根据回归系数和特征向量来计算 Sigmoid 的值，大于0.5函数返回1，否则返回0
    Args:
        inX -- 特征向量，features
        weights -- 根据梯度下降/随机梯度下降 计算得到的回归系数
    Returns:
        如果 prob 计算大于 0.5 函数返回 1
        否则返回 0
    '''
    prob = sigmoid(sum(inX * weights))
    if prob > 0.5: return 1.0
    else: return 0.0

# 打开测试集和训练集,并对数据进行格式化处理
def colicTest():
    '''
    Desc:
        打开测试集和训练集，并对数据进行格式化处理
    Args:
        None
    Returns:
        errorRate -- 分类错误率
    '''
    frTrain = open('input/5.Logistic/horseColicTraining.txt')
    frTest = open('input/5.Logistic/horseColicTest.txt')
    trainingSet = []
    trainingLabels = []
    # 解析训练数据集中的数据特征和Labels
    # trainingSet 中存储训练数据集的特征，trainingLabels 存储训练数据集的样本对应的分类标签
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    # 使用 改进后的 随机梯度下降算法 求得在此数据集上的最佳回归系数 trainWeights
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 500)
    # trainWeights = stocGradAscent0(array(trainingSet), trainingLabels)
    errorCount = 0
    numTestVec = 0.0
    # 读取 测试数据集 进行测试，计算分类错误的样本条数和最终的错误率
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(array(lineArr), trainWeights)) != int(
                currLine[21]):
            errorCount += 1
    errorRate = (float(errorCount) / numTestVec)
    print "the error rate of this test is: %f" % errorRate
    return errorRate

# 调用 colicTest() 10次并求结果的平均值
def multiTest():
    numTests = 10
    errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print "after %d iterations the average error rate is: %f" % (numTests, errorSum / float(numTests))

if __name__ == "__main__":
    simpleTest()
    # multiTest()

reason: from numpy