text:
#
# Copyright (C) 2020 Pico Technology Ltd. See LICENSE file for terms.
#
# PS6000 A BLOCK MODE EXAMPLE
# This example opens a 6000a driver device, sets up two channels and a trigger then collects a block of data.
# This data is then plotted as mV against time in ns.

import ctypes
import numpy as np
from picosdk.ps6000a import ps6000a as ps
from picosdk.PicoDeviceEnums import picoEnum as enums
import matplotlib.pyplot as plt
from picosdk.functions import adc2mV, assert_pico_ok

# Create chandle and status ready for use
chandle = ctypes.c_int16()
status = {}

# Open 6000 A series PicoScope
# returns handle to chandle for use in future API functions
resolution = enums.PICO_DEVICE_RESOLUTION["PICO_DR_8BIT"]
status["openunit"] = ps.ps6000aOpenUnit(ctypes.byref(chandle), None, resolution)
assert_pico_ok(status["openunit"])

# Set channel A on
# handle = chandle
channelA = enums.PICO_CHANNEL["PICO_CHANNEL_A"]
coupling = enums.PICO_COUPLING["PICO_DC"]
channelRange = 7
# analogueOffset = 0 V
bandwidth = enums.PICO_BANDWIDTH_LIMITER["PICO_BW_FULL"]
status["setChannelA"] = ps.ps6000aSetChannelOn(chandle, channelA, coupling, channelRange, 0, bandwidth)
assert_pico_ok(status["setChannelA"])

# set channel B-H off
for x in range (1, 7, 1):
    channel = x
    status["setChannel",x] = ps.ps6000aSetChannelOff(chandle,channel)
    assert_pico_ok(status["setChannel",x])

# Set simple trigger on channel A, 1 V rising with 1 s autotrigger
# handle = chandle
# enable = 1
source = channelA
# threshold = 1000 mV
direction = enums.PICO_THRESHOLD_DIRECTION["PICO_RISING"]
# delay = 0 s
# autoTriggerMicroSeconds = 1000000 us
status["setSimpleTrigger"] = ps.ps6000aSetSimpleTrigger(chandle, 1, source, 1000, direction, 0, 1000000)
assert_pico_ok(status["setSimpleTrigger"])

# Get fastest available timebase
# handle = chandle
enabledChannelFlags = enums.PICO_CHANNEL_FLAGS["PICO_CHANNEL_A_FLAGS"]
timebase = ctypes.c_uint32(0)
timeInterval = ctypes.c_double(0)
# resolution = resolution
status["getMinimumTimebaseStateless"] = ps.ps6000aGetMinimumTimebaseStateless(chandle, enabledChannelFlags, ctypes.byref(timebase), ctypes.byref(timeInterval), resolution)
print("timebase = ", timebase.value)
print("sample interval =", timeInterval.value, "s")

# Set number of samples to be collected
noOfPreTriggerSamples = 500000
noOfPostTriggerSamples = 1000000
nSamples = noOfPostTriggerSamples + noOfPreTriggerSamples

# Create buffers
bufferAMax = (ctypes.c_int16 * nSamples)()
bufferAMin = (ctypes.c_int16 * nSamples)() # used for downsampling which isn't in the scope of this example

# Set data buffers
# handle = chandle
# channel = channelA
# bufferMax = bufferAMax
# bufferMin = bufferAMin
# nSamples = nSamples
dataType = enums.PICO_DATA_TYPE["PICO_INT16_T"]
waveform = 0
downSampleMode = enums.PICO_RATIO_MODE["PICO_RATIO_MODE_RAW"]
clear = enums.PICO_ACTION["PICO_CLEAR_ALL"]
add = enums.PICO_ACTION["PICO_ADD"]
action = clear|add # PICO_ACTION["PICO_CLEAR_WAVEFORM_CLEAR_ALL"] | PICO_ACTION["PICO_ADD"]
status["setDataBuffers"] = ps.ps6000aSetDataBuffers(chandle, channelA, ctypes.byref(bufferAMax), ctypes.byref(bufferAMin), nSamples, dataType, waveform, downSampleMode, action)
assert_pico_ok(status["setDataBuffers"])

# Run block capture
# handle = chandle
# timebase = timebase
timeIndisposedMs = ctypes.c_double(0)
# segmentIndex = 0
# lpReady = None   Using IsReady rather than a callback
# pParameter = None
status["runBlock"] = ps.ps6000aRunBlock(chandle, noOfPreTriggerSamples, noOfPostTriggerSamples, timebase, ctypes.byref(timeIndisposedMs), 0, None, None)
assert_pico_ok(status["runBlock"])

# Check for data collection to finish using ps6000aIsReady
ready = ctypes.c_int16(0)
check = ctypes.c_int16(0)
while ready.value == check.value:
    status["isReady"] = ps.ps6000aIsReady(chandle, ctypes.byref(ready))

# Get data from scope
# handle = chandle
# startIndex = 0
noOfSamples = ctypes.c_uint64(nSamples)
# downSampleRatio = 1
# segmentIndex = 0
overflow = ctypes.c_int16(0)
status["getValues"] = ps.ps6000aGetValues(chandle, 0, ctypes.byref(noOfSamples), 1, downSampleMode, 0, ctypes.byref(overflow))
assert_pico_ok(status["getValues"])

# get max ADC value
# handle = chandle
minADC = ctypes.c_int16()
maxADC = ctypes.c_int16()
status["getAdcLimits"] = ps.ps6000aGetAdcLimits(chandle, resolution, ctypes.byref(minADC), ctypes.byref(maxADC))
assert_pico_ok(status["getAdcLimits"])

# convert ADC counts data to mV
adc2mVChAMax =  adc2mV(bufferAMax, channelRange, maxADC)

# Create time data
time = np.linspace(0, (nSamples -1) * timeInterval.value * 1000000000, nSamples)

# plot data from channel A and B
plt.plot(time, adc2mVChAMax[:])
plt.xlabel('Time (ns)')
plt.ylabel('Voltage (mV)')
plt.show()

# Close the scope
status["closeunit"] = ps.ps6000aCloseUnit(chandle)
assert_pico_ok(status["closeunit"])

print(status)

max_stars_repo_name: LauritzRaisch/picosdk-python-wrappers

max_issues_repo_name: LauritzRaisch/picosdk-python-wrappers

max_forks_repo_name: LauritzRaisch/picosdk-python-wrappers


text:
# Copyright 2021 Sony Corporation.
# Copyright 2021 Sony Group Corporation.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import numpy as np

import nnabla as nn
import nnabla.functions as F

from model import PretrainedVgg19

def perceptual_loss(pyramide_real, pyramide_fake, scales, weights, vgg_param_path):
    """
        Compute Perceptual Loss using VGG19 as a feature extractor.
    """
    vgg19 = PretrainedVgg19(param_path=vgg_param_path)
    variable_not_exist = True
    for scale in scales:
        x_vgg = vgg19(pyramide_fake[f'prediction_{scale}'])
        y_vgg = vgg19(pyramide_real[f'prediction_{scale}'])

        for i, weight in enumerate(weights):
            value = F.mean(F.absolute_error(x_vgg[i], y_vgg[i]))
            if variable_not_exist:
                loss = weight * value
                variable_not_exist = False
            else:
                loss += weight * value
    return loss

def lsgan_loss(real, weight, fake=None):
    if fake:
        loss = weight * F.mean(F.squared_error(F.constant(1, real.shape), real)
                               + F.pow_scalar(fake, 2))
    else:
        loss = weight * \
            F.mean(F.squared_error(F.constant(1, real.shape), real))
    return loss

def feature_matching_loss(dis_maps_real, dis_maps_fake, model_params, weights):
    variable_not_exist = True
    for j, scale in enumerate(model_params['discriminator_params']['scales']):
        key = f"feature_maps_{scale}".replace('.', '-')
        for i, (a, b) in enumerate(zip(dis_maps_real[key], dis_maps_fake[key])):
            if weights[i] == 0:
                continue
            if variable_not_exist:
                loss = F.mean(F.absolute_error(a, b)) * weights[i]
                variable_not_exist = False
            else:
                loss += F.mean(F.absolute_error(a, b)) * weights[i]
    return loss

def equivariance_value_loss(kp_driving_value, warped_kp_value, weight):
    value_loss = F.mean(F.absolute_error(kp_driving_value, warped_kp_value))
    loss = weight * value_loss
    return loss

def equivariance_jacobian_loss(kp_driving_jacobian,
                               arithmetic_jacobian,
                               trans_kp_jacobian, weight):
    jacobian_transformed = F.batch_matmul(arithmetic_jacobian,
                                          trans_kp_jacobian)

    normed_driving = F.reshape(
                        F.batch_inv(
                            F.reshape(kp_driving_jacobian,
                                      (-1,) + kp_driving_jacobian.shape[-2:])),
                        kp_driving_jacobian.shape)

    normed_transformed = jacobian_transformed
    value = F.batch_matmul(normed_driving, normed_transformed)

    eye = nn.Variable.from_numpy_array(np.reshape(np.eye(2), (1, 1, 2, 2)))

    jacobian_loss = F.mean(F.absolute_error(eye, value))
    loss = weight * jacobian_loss
    return loss

ext: py

max_stars_repo_name: AaratiAkkapeddi/nnabla-examples

max_issues_repo_name: AaratiAkkapeddi/nnabla-examples

max_forks_repo_name: AaratiAkkapeddi/nnabla-examples

reason: import numpy


text:
#!/usr/bin/env python

"""
module for manipulating netcdf particle files

This is a test case for working with what hopefully will be a CF standard

"""
from __future__ import division #Change the / operator to ensure true division throughout (Zelenke).
from datetime import datetime

import numpy as np

import netCDF4

class particle_trajectory:
    def __init__ (self):

        # lots of defaults:
        self.Trajectory = []
        self.timesteps = []
        self.time_units = "hours since 2010-01-01 00:00:00"
        self.attributes = {'Conventions' : "CF-1.6",
                           'title' : "Sample data/file for particle trajectory format",
                           'institution' : "NOAA Emergency Response Division",
                           'source' : "Example Data",
                           'history' : "Evolved with discussion on CF-metadata listserve",
                           'references' : '',
                           'comment' : 'Some simple test data',
                           'CF:featureType' : "particle_trajectory" ,
                           }

#    dimensions:
#   record = UNLIMITED ; // (15768 currently)
#   time = 145 ;
#variables:
#
#       int number_records_per_timestep(time) ;
#       number_records_per_timestep:units = "1" ;
#       number_records_per_timestep:long_name = "number spillets per current timestep" ;
#       number_records_per_timestep:CF\:ragged_row_count = "record" ;
#   double time(record) ;
#       time:long_name = "forecast time after simulation start" ;
#       time:standard_name = "forecast_reference_time" ;
#       time:units = "days since 2010-05-25" ;
#       time:calendar = "gregorian" ;
#   float latitude(record) ;
#       latitude:long_name = "contaminant slick latitude" ;
#       latitude:standard_name = "latitude" ;
#       latitude:units = "degrees_north" ;
#   float longitude(record) ;

    def write(self, filename):
        nc = netCDF4.Dataset(filename,'w',format='NETCDF3_CLASSIC')

        # Global Attributes
        for attr, val in self.attributes.items():
            nc.setncattr(attr, val)
            setattr(nc, 'creation_date', datetime.now().isoformat())

        # add Dimensions
        print "num timesteps:", len(self.Trajectory)
        nc.createDimension('time', len(self.Trajectory))
        nc.createDimension('data', None)

        # create variables
        print "creating variables"
        ## fixme: should be able to create these from Trajectory dtype.
        Variables = { 'time' : nc.createVariable('time',  'f4',('time',)),
                      'particle_count': nc.createVariable('particle_count', 'i4',('time')),
                      'longitude'  : nc.createVariable('longitude',   'f4',('data',)),
                      'latitude'  : nc.createVariable('latitude',   'f4',('data',)),
                      'depth': nc.createVariable('depth', 'f4',('data',)),
                      'mass' : nc.createVariable('mass',  'f4',('data',)),
                      'flag' : nc.createVariable('flag',  'i1',('data',)),
                      'id'   : nc.createVariable('id'  ,  'i4',('data',)),
                      }
        #diameter  = nc.createVariable('diameter', 'f4',('time','particle'))
        #rise_velocity  = nc.createVariable('rise_velocity', 'f4',('time','particle'))
        #vol = nc.createVariable('vol', 'f4',('time','particle'))
        #density  = nc.createVariable('density', 'f4',('time','particle'))
        #release_time  = nc.createVariable('release_time', 'f4',('time','particle'))

        # add data
        print "adding data"
        ##fixme -- should be able to add on by time -- how?
        for field in self.Trajectory[0].dtype.names:
            data = []
            for t in range(len(self.Trajectory)):
                data.extend(self.Trajectory[t][field])
                Variables['particle_count'][t] = len(self.Trajectory[t])
            Variables[field][:] = data

        # time
        Variables['time'][:] = netCDF4.date2num(self.timesteps, self.time_units)

        ## add attributes:
        Variable_attributes = {}
        Variable_attributes['time']  = {'long_name' : "Time",
                                        'standard_name' : "time",
                                        'units' : self.time_units
                                        }

        Variable_attributes['particle_count'] = {'units': "1",
                                            'long_name': "number particles in given timestep",
                                            'CF:ragged_row_count' : "data",
                                            }

        Variable_attributes['lon']  = {'long_name' : "Longitude of the particle",
                                       }
#            lon.standard_name = "longitude_particle"
#            lon.Fill_Value = -99999.
#            lon.missing_value = -99999.
#            lon.units = "degrees"

        Variable_attributes['lat']  = {'long_name' : "Latitude of the particle",
                                       }

#            lat.long_name = "Latitude of the particle"
#            lat.standard_name = "latitude_particle"
#            lat.Fill_Value = -99999.
#            lat.missing_value = -99999.
#            lat.units = "degrees"

#            depth.long_name = "Depth of the particles"
#            depth.standard_name = "depth_particle"
#            depth.Fill_Value = -99999.
#            depth.missing_value = -99999.
#            depth.units = "meter"
        for var, attrs in Variable_attributes.items():
            for att, val in attrs.items():
                Variables[var].setncattr(att, val)
        nc.close()

#DSM: Wrapping of particle data generated by GNOME simulation
class nc_particle_file:
    """
    class to wrap a netcdf particle file
    """

    def __init__(self, nc):

        self.nc = nc

        time = nc.variables['time']
        units = time.getncattr('units')
        self.times = netCDF4.num2date(time, units)
        self.time_units = units

        #Defined mass in the same way as done above for time (Zelenke).
        mass = nc.variables['mass']
        units_mass = mass.getncattr('units')
        self.mass = mass
        self.mass_units = units_mass

        #print nc.variables

        self.particle_count = nc.variables['particle_count']
        # build the index:
        self.data_index = np.zeros((len(self.times)+1,), dtype=np.int32 )
        self.data_index[1:] = np.cumsum(self.particle_count)

        #print self.times
        #print self.particle_count
        #print self.data_index

    def get_all_timesteps(self, variables=['latitude','longitude']):
         """
         returns the requested variables data from a given timestep as a
         dictionary keyed by the variable names
         """
         data = {}
         for var in variables:
             data[var] = self.nc.variables[var][:]
         return data

    def get_timestep(self, timestep, variables=['latitude','longitude']):
        """
        returns the requested variables data from a given timestep as a
        dictionary keyed by the variable names
        """
        data = {}
        for var in variables:
            ind1 = self.data_index[timestep]
            ind2 = self.data_index[timestep+1]
            data[var] = self.nc.variables[var][ind1:ind2]
        return data

    def get_individual_trajectory(self, particle_id, vars=['latitude','longitude']):
        """
        returns the requested variables from trajectory of an individual particle

        note: this is very inefficient -- it has to read the entire file to get it.
        """
        indexes = np.where(self.nc.variables['id'] == particle_id)
        #print self.nc.variables['id'][:]
        #print indexes

max_stars_repo_name: dylanrighi/GnomeTools

max_issues_repo_name: dylanrighi/GnomeTools

max_forks_repo_name: dylanrighi/GnomeTools


text:
import pandas as pd
import numpy as np
np.random.seed(99)
from sklearn.model_selection import train_test_split
from sklearn.model_selection import KFold
from sklearn.model_selection import GridSearchCV
from sklearn.multioutput import MultiOutputClassifier, MultiOutputRegressor
from sklearn.multiclass import OneVsRestClassifier
import xgboost as xgb
from xgboost.sklearn import XGBClassifier
from xgboost.sklearn import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.multiclass import OneVsRestClassifier
from sklearn.preprocessing import LabelEncoder
import lightgbm as lgbm
from sklearn.model_selection import KFold, cross_val_score,StratifiedKFold
import seaborn as sns
from sklearn.preprocessing import OneHotEncoder, LabelEncoder, label_binarize
import csv
import re
from xgboost import XGBRegressor, XGBClassifier
from sklearn.metrics import mean_squared_log_error, mean_squared_error,balanced_accuracy_score
from scipy import stats
from sklearn.model_selection import RandomizedSearchCV
import scipy as sp
import time
import copy
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from collections import Counter, defaultdict
import pdb
#################  All these imports are needed for the pipeline #######
import time
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import OneHotEncoder
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.linear_model import LogisticRegression
from sklearn.compose import make_column_transformer
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import LabelEncoder, LabelBinarizer
from sklearn.base import BaseEstimator, TransformerMixin #gives fit_transform method for free
import pdb
from sklearn.base import TransformerMixin
from collections import defaultdict
from sklearn.preprocessing import MaxAbsScaler
from sklearn.preprocessing import FunctionTransformer
###################################################################################################
#############   This is where you import from other Auto_TS modules ############
from ..utils import My_LabelEncoder, My_LabelEncoder_Pipe
from ..utils import left_subtract
#################################################################################
def complex_XGBoost_model(X_train, y_train, X_test, log_y=False, GPU_flag=False,
                                scaler = '', enc_method='label', n_splits=5,
                                num_boost_round=1000, verbose=-1):
    """
    This model is called complex because it handle multi-label, mulit-class datasets which XGBoost ordinarily cant.
    Just send in X_train, y_train and what you want to predict, X_test
    It will automatically split X_train into multiple folds (10) and train and predict each time on X_test.
    It will then use average (or use mode) to combine the results and give you a y_test.
    It will automatically detect modeltype as "Regression" or 'Classification'
    It will also add MultiOutputClassifier and MultiOutputRegressor to multi_label problems.
    The underlying estimators in all cases is XGB. So you get the best of both worlds.

    Inputs:
    ------------
    X_train: pandas dataframe only: do not send in numpy arrays. This is the X_train of your dataset.
    y_train: pandas Series or DataFrame only: do not send in numpy arrays. This is the y_train of your dataset.
    X_test: pandas dataframe only: do not send in numpy arrays. This is the X_test of your dataset.
    log_y: default = False: If True, it means use the log of the target variable "y" to train and test.
    GPU_flag: if your machine has a GPU set this flag and it will use XGBoost GPU to speed up processing.
    scaler : default is StandardScaler(). But you can send in MinMaxScaler() as input to change it or any other scaler.
    enc_method: default is 'label' encoding. But you can choose 'glmm' as an alternative. But those are the only two.
    verbose: default = 0. Choosing 1 will give you lot more output.

    Outputs:
    ------------
    y_preds: Predicted values for your X_XGB_test dataframe.
        It has been averaged after repeatedly predicting on X_XGB_test. So likely to be better than one model.
    """
    X_XGB = copy.deepcopy(X_train)
    Y_XGB = copy.deepcopy(y_train)
    X_XGB_test = copy.deepcopy(X_test)
    ####################################
    start_time = time.time()
    top_num = 10
    if isinstance(Y_XGB, pd.Series):
        targets = [Y_XGB.name]
    else:
        targets = Y_XGB.columns.tolist()
    if len(targets) == 1:
        multi_label = False
        if isinstance(Y_XGB, pd.DataFrame):
            Y_XGB = pd.Series(Y_XGB.values.ravel(),name=targets[0], index=Y_XGB.index)
    else:
        multi_label = True
    modeltype, _ = analyze_problem_type(Y_XGB, targets)
    columns =  X_XGB.columns
    ##### Now continue with scaler pre-processing ###########
    if isinstance(scaler, str):
        if not scaler == '':
            scaler = scaler.lower()
    if scaler == 'standard':
        scaler = StandardScaler()
    elif scaler == 'minmax':
        scaler = MinMaxScaler()
    else:
        scaler = StandardScaler()
    #########     G P U     P R O C E S S I N G      B E G I N S    ############
    ###### This is where we set the CPU and GPU parameters for XGBoost
    if GPU_flag:
        GPU_exists = check_if_GPU_exists()
    else:
        GPU_exists = False
    #####   Set the Scoring Parameters here based on each model and preferences of user ###
    cpu_params = {}
    param = {}
    cpu_params['tree_method'] = 'hist'
    cpu_params['gpu_id'] = 0
    cpu_params['updater'] = 'grow_colmaker'
    cpu_params['predictor'] = 'cpu_predictor'
    if GPU_exists:
        param['tree_method'] = 'gpu_hist'
        param['gpu_id'] = 0
        param['updater'] = 'grow_gpu_hist' #'prune'
        param['predictor'] = 'gpu_predictor'
        print('    Hyper Param Tuning XGBoost with GPU parameters. This will take time. Please be patient...')
    else:
        param = copy.deepcopy(cpu_params)
        print('    Hyper Param Tuning XGBoost with CPU parameters. This will take time. Please be patient...')
    #################################################################################
    if modeltype == 'Regression':
        if log_y:
            Y_XGB.loc[Y_XGB==0] = 1e-15  ### just set something that is zero to a very small number

    #########  Now set the number of rows we need to tune hyper params ###
    scoreFunction = { "precision": "precision_weighted","recall": "recall_weighted"}
    random_search_flag =  True

      #### We need a small validation data set for hyper-param tuning #########################
    hyper_frac = 0.2
    #### now select a random sample from X_XGB ##

    if modeltype == 'Regression':
        X_train, X_valid, Y_train, Y_valid = train_test_split(X_XGB, Y_XGB, test_size=hyper_frac,
                            random_state=999)
    else:
        X_train, X_valid, Y_train, Y_valid = train_test_split(X_XGB, Y_XGB, test_size=hyper_frac,
                            random_state=999, stratify = Y_XGB)

    ######  This step is needed for making sure y is transformed to log_y ####################
    if modeltype == 'Regression' and log_y:
            Y_train = np.log(Y_train)
            Y_valid = np.log(Y_valid)

    #### First convert test data into numeric using train data ###
    X_train, Y_train, X_valid, Y_valid, scaler = data_transform(X_train, Y_train, X_valid, Y_valid,
                                modeltype, multi_label, scaler=scaler, enc_method=enc_method)

    ######  Time to hyper-param tune model using randomizedsearchcv and partial train data #########
    num_boost_round = xgbm_model_fit(random_search_flag, X_train, Y_train, X_valid, Y_valid, modeltype,
                         multi_label, log_y, num_boost_round=num_boost_round)

    #### First convert test data into numeric using train data ###############################
    if not isinstance(X_XGB_test, str):
        x_train, y_train, x_test, _, _ = data_transform(X_XGB, Y_XGB, X_XGB_test, "",
                                modeltype, multi_label, scaler=scaler, enc_method=enc_method)

    ######  Time to train the hyper-tuned model on full train data ##########################
    random_search_flag = False
    model = xgbm_model_fit(random_search_flag, x_train, y_train, x_test, "", modeltype,
                                multi_label, log_y, num_boost_round=num_boost_round)

    #############  Time to get feature importances based on full train data   ################

    if multi_label:
        for i,target_name in enumerate(targets):
            each_model = model.estimators_[i]
            imp_feats = dict(zip(x_train.columns, each_model.feature_importances_))
            importances = pd.Series(imp_feats).sort_values(ascending=False)[:top_num].values
            important_features = pd.Series(imp_feats).sort_values(ascending=False)[:top_num].index.tolist()
            print('Top 10 features for {}: {}'.format(target_name, important_features))
    else:
        imp_feats = model.get_score(fmap='', importance_type='gain')
        importances = pd.Series(imp_feats).sort_values(ascending=False)[:top_num].values
        important_features = pd.Series(imp_feats).sort_values(ascending=False)[:top_num].index.tolist()
        print('Top 10 features:\n%s' %important_features[:top_num])
        #######  order this in the same order in which they were collected ######
        feature_importances = pd.DataFrame(importances,
                                           index = important_features,
                                            columns=['importance'])

    ######  Time to consolidate the predictions on test data ################################
    if not multi_label and not isinstance(X_XGB_test, str):
        x_test = xgb.DMatrix(x_test)
    if isinstance(X_XGB_test, str):
        print('No predictions since X_XGB_test is empty string. Returning...')
        return {}
    if modeltype == 'Regression':
        if not isinstance(X_XGB_test, str):
            if log_y:
                pred_xgbs = np.exp(model.predict(x_test))
            else:
                pred_xgbs = model.predict(x_test)
            #### if there is no test data just return empty strings ###
        else:
            pred_xgbs = []
    else:
        if multi_label:
            pred_xgbs = model.predict(x_test)
            pred_probas = model.predict_proba(x_test)
        else:
            pred_probas = model.predict(x_test)
            if modeltype =='Multi_Classification':
                pred_xgbs = pred_probas.argmax(axis=1)
            else:
                pred_xgbs = (pred_probas>0.5).astype(int)
    ##### once the entire model is trained on full train data ##################
    print('    Time taken for training XGBoost on entire train data (in minutes) = %0.1f' %(
             (time.time()-start_time)/60))
    if multi_label:
        for i,target_name in enumerate(targets):
            each_model = model.estimators_[i]
            xgb.plot_importance(each_model, importance_type='gain', title='XGBoost model feature importances for %s' %target_name)
    else:
        xgb.plot_importance(model, importance_type='gain', title='XGBoost final model feature importances')
    print('Returning the following:')
    print('    Model = %s' %model)
    print('    Scaler = %s' %scaler)
    if modeltype == 'Regression':
        if not isinstance(X_XGB_test, str):
            print('    (3) sample predictions:%s' %pred_xgbs[:3])
        return (pred_xgbs, scaler, model)
    else:
        if not isinstance(X_XGB_test, str):
            print('    (3) sample predictions (may need to be transformed to original labels):%s' %pred_xgbs[:3])
            print('    (3) sample predicted probabilities:%s' %pred_probas[:3])
        return (pred_xgbs, scaler, model)
##############################################################################################
import xgboost as xgb
def xgbm_model_fit(random_search_flag, x_train, y_train, x_test, y_test, modeltype,
                         multi_label, log_y, num_boost_round=100):
    start_time = time.time()
    if multi_label and not random_search_flag:
        model = num_boost_round
    else:
        rand_params = {
            'learning_rate': sp.stats.uniform(scale=1),
            'gamma': sp.stats.randint(0, 100),
            'n_estimators': sp.stats.randint(100,500),
            "max_depth": sp.stats.randint(3, 15),
                }

    if modeltype == 'Regression':
        objective = 'reg:squarederror'
        eval_metric = 'rmse'
        shuffle = False
        stratified = False
        num_class = 0
        score_name = 'Score'
        scale_pos_weight = 1
    else:
        if modeltype =='Binary_Classification':
            objective='binary:logistic'
            eval_metric = 'error' ## dont foolishly change to auc or aucpr since it doesnt work in finding feature imps later
            shuffle = True
            stratified = True
            num_class = 1
            score_name = 'Error Rate'
            scale_pos_weight = get_scale_pos_weight(y_train)
        else:
            objective = 'multi:softprob'
            eval_metric = 'merror'  ## dont foolishly change to auc or aucpr since it doesnt work in finding feature imps later
            shuffle = True
            stratified = True
            if multi_label:
                num_class = y_train.nunique().max()
            else:
                if isinstance(y_train, np.ndarray):
                    num_class = np.unique(y_train).max() + 1
                elif isinstance(y_train, pd.Series):
                    num_class = y_train.nunique()
                else:
                    num_class = y_train.nunique().max()
            score_name = 'Multiclass Error Rate'
            scale_pos_weight = 1  ### use sample_weights in multi-class settings ##
    ######################################################
    final_params = {
          'booster' :'gbtree',
          'colsample_bytree': 0.5,
          'alpha': 0.015,
          'gamma': 4,
          'learning_rate': 0.01,
          'max_depth': 8,
          'min_child_weight': 2,
          'reg_lambda': 0.5,
          'subsample': 0.7,
          'random_state': 99,
          'objective': objective,
          'eval_metric': eval_metric,
          'verbosity': 0,
          'n_jobs': -1,
          'scale_pos_weight':scale_pos_weight,
          'num_class': num_class,
          'silent': True
            }
    #######  This is where we split into single and multi label ############
    if multi_label:
        ######   This is for Multi_Label problems ############
        rand_params = {'estimator__learning_rate':[0.1, 0.5, 0.01, 0.05],
          'estimator__n_estimators':[50, 100, 150, 200, 250],
          'estimator__gamma':[2, 4, 8, 16, 32],
          'estimator__max_depth':[3, 5, 8, 12],
          }
        if random_search_flag:
            if modeltype == 'Regression':
                clf = XGBRegressor(n_jobs=-1, random_state=999, max_depth=6)
                clf.set_params(**final_params)
                model = MultiOutputRegressor(clf, n_jobs=-1)
            else:
                clf = XGBClassifier(n_jobs=-1, random_state=999, max_depth=6)
                clf.set_params(**final_params)
                model = MultiOutputClassifier(clf, n_jobs=-1)
            if modeltype == 'Regression':
                scoring = 'neg_mean_squared_error'
            else:
                scoring = 'precision'
            model = RandomizedSearchCV(model,
                       param_distributions = rand_params,
                       n_iter = 15,
                       return_train_score = True,
                       random_state = 99,
                       n_jobs=-1,
                       cv = 3,
                       refit=True,
                       scoring = scoring,
                       verbose = False)
            model.fit(x_train, y_train)
            print('Time taken for Hyper Param tuning of multi_label XGBoost (in minutes) = %0.1f' %(
                                            (time.time()-start_time)/60))
            cv_results = pd.DataFrame(model.cv_results_)
            print('Mean cross-validated test %s = %0.04f' %(score_name, cv_results['mean_test_score'].mean()))
            ### In this case, there is no boost rounds so just return the default num_boost_round
            return model.best_estimator_
        else:
            try:
                model.fit(x_train, y_train)
            except:
                print('Multi_label XGBoost model is crashing during training. Please check your inputs and try again...')
            return model
    else:
        #### This is for Single Label Problems #############
        if modeltype == 'Multi_Classification':
            wt_array = get_sample_weight_array(y_train)
            dtrain = xgb.DMatrix(x_train, label=y_train, weight=wt_array)
        else:
            dtrain = xgb.DMatrix(x_train, label=y_train)
        ########   Now let's perform randomized search to find best hyper parameters ######
        if random_search_flag:
            cv_results = xgb.cv(final_params, dtrain, num_boost_round=num_boost_round, nfold=5,
                stratified=stratified, metrics=eval_metric, early_stopping_rounds=10, seed=999, shuffle=shuffle)
            # Update best eval_metric
            best_eval = 'test-'+eval_metric+'-mean'
            mean_mae = cv_results[best_eval].min()
            boost_rounds = cv_results[best_eval].argmin()
            print("Cross-validated %s = %0.3f in num rounds = %s" %(score_name, mean_mae, boost_rounds))
            print('Time taken for Hyper Param tuning of XGBoost (in minutes) = %0.1f' %(
                                                (time.time()-start_time)/60))
            return boost_rounds
        else:
            try:
                model = xgb.train(
                    final_params,
                    dtrain,
                    num_boost_round=num_boost_round,
                    verbose_eval=False,
                )
            except:
                print('XGBoost model is crashing. Please check your inputs and try again...')
            return model
####################################################################################
# Calculate class weight
from sklearn.utils.class_weight import compute_class_weight
import copy
from collections import Counter
def find_rare_class(classes, verbose=0):
    ######### Print the % count of each class in a Target variable  #####
    """
    Works on Multi Class too. Prints class percentages count of target variable.
    It returns the name of the Rare class (the one with the minimum class member count).
    This can also be helpful in using it as pos_label in Binary and Multi Class problems.
    """
    counts = OrderedDict(Counter(classes))
    total = sum(counts.values())
    if verbose >= 1:
        print('       Class  -> Counts -> Percent')
        sorted_keys = sorted(counts.keys())
        for cls in sorted_keys:
            print("%12s: % 7d  ->  % 5.1f%%" % (cls, counts[cls], counts[cls]/total*100))
    if type(pd.Series(counts).idxmin())==str:
        return pd.Series(counts).idxmin()
    else:
        return int(pd.Series(counts).idxmin())
###################################################################################
def get_sample_weight_array(y_train):
    y_train = copy.deepcopy(y_train)
    if isinstance(y_train, np.ndarray):
        y_train = pd.Series(y_train)
    elif isinstance(y_train, pd.Series):
        pass
    elif isinstance(y_train, pd.DataFrame):
        ### if it is a dataframe, return only if it s one column dataframe ##
        y_train = y_train.iloc[:,0]
    else:
        ### if you cannot detect the type or if it is a multi-column dataframe, ignore it
        return None
    classes = np.unique(y_train)
    class_weights = compute_class_weight('balanced', classes=classes, y=y_train)
    if len(class_weights[(class_weights < 1)]) > 0:
        ### if the weights are less than 1, then divide them until the lowest weight is 1.
        class_weights = class_weights/min(class_weights)
    else:
        class_weights = (class_weights)
    ### even after you change weights if they are all below 1.5 do this ##
    #if (class_weights<=1.5).all():
    #    class_weights = np.around(class_weights+0.49)
    class_weights = class_weights.astype(int)
    wt = dict(zip(classes, class_weights))

    ### Map class weights to corresponding target class values
    ### You have to make sure class labels have range (0, n_classes-1)
    wt_array = y_train.map(wt)
    #set(zip(y_train, wt_array))

    # Convert wt series to wt array
    wt_array = wt_array.values
    return wt_array
###############################################################################
from collections import OrderedDict
def get_scale_pos_weight(y_input):
    y_input = copy.deepcopy(y_input)
    if isinstance(y_input, np.ndarray):
        y_input = pd.Series(y_input)
    elif isinstance(y_input, pd.Series):
        pass
    elif isinstance(y_input, pd.DataFrame):
        ### if it is a dataframe, return only if it s one column dataframe ##
        y_input = y_input.iloc[:,0]
    else:
        ### if you cannot detect the type or if it is a multi-column dataframe, ignore it
        return None
    classes = np.unique(y_input)
    rare_class = find_rare_class(y_input)
    xp = Counter(y_input)
    class_weights = compute_class_weight('balanced', classes=classes, y=y_input)

    if len(class_weights[(class_weights < 1)]) > 0:
        ### if the weights are less than 1, then divide them until the lowest weight is 1.
        class_weights = class_weights/min(class_weights)
    else:
        class_weights = (class_weights)
    ### even after you change weights if they are all below 1.5 do this ##
    #if (class_weights<=1.5).all():
    #    class_weights = np.around(class_weights+0.49)

    class_weights = class_weights.astype(int)
    class_weights[(class_weights<1)]=1
    class_rows = class_weights*[xp[x] for x in classes]
    class_rows = class_rows.astype(int)
    class_weighted_rows = dict(zip(classes,class_weights))
    rare_class_weight = class_weighted_rows[rare_class]
    print('    For class %s, weight = %s' %(rare_class, rare_class_weight))
    return rare_class_weight
#########################################################################################################
###########################################################################################
from collections import defaultdict
from collections import OrderedDict
from sklearn.impute import SimpleImputer
def data_transform(X_train, Y_train, X_test="", Y_test="", modeltype='Classification',
            multi_label=False, enc_method='label', scaler=""):

    #### All these are needed for transforming cat variables and building a pipeline ###
    imp_constant = SimpleImputer(strategy='constant', fill_value='missing')
    ohe = OneHotEncoder()
    imp_ohe = make_pipeline(imp_constant, ohe)
    vect = CountVectorizer()
    imp = SimpleImputer()
    le = My_LabelEncoder()

    def drop_second_col(Xt):
        ### This deletes the 2nd column. Hence col number=1 and axis=1 ###
        return np.delete(Xt, 1, 1)

    ####  This is where we define the Pipeline for cat encoders and Label Encoder ############
    lep = My_LabelEncoder_Pipe()
    drop_second_col_func = FunctionTransformer(drop_second_col)
    ### lep_one uses My_LabelEncoder to first label encode and then drop the second unused column ##
    lep_one = make_pipeline(lep, drop_second_col_func)

    ### if you drop remainder variables, then leftovervars is not needed.
    ### If you passthrough remainder variables, then leftovers must be included
    remainder = 'drop'

    ### If you choose MaxAbsScaler, then NaNs which were Label Encoded as -1 are preserved as - (negatives). This is fantastic.
    ### If you choose StandardScaler or MinMaxScaler, the integer values become stretched as if they are far
    ###    apart when in reality they are close. So avoid it for now.
    scaler = MaxAbsScaler()
    #scaler = StandardScaler()
    ##### First make sure that the originals are not modified ##########
    X_train_encoded = copy.deepcopy(X_train)
    X_test_encoded = copy.deepcopy(X_test)
    ##### Use My_Label_Encoder to transform label targets if needed #####
    if multi_label:
        if modeltype != 'Regression':
            targets = Y_train.columns
            Y_train_encoded = copy.deepcopy(Y_train)
            for each_target in targets:
                mlb = My_LabelEncoder()
                if not isinstance(Y_train, str):
                    Y_train_encoded[each_target] = mlb.fit_transform(Y_train[each_target])
                else:
                    Y_train_encoded = copy.deepcopy(Y_train)
                if not isinstance(Y_test, str):
                    Y_test_encoded= mlb.transform(Y_test)
                else:
                    Y_test_encoded = copy.deepcopy(Y_test)
        else:
            Y_train_encoded = copy.deepcopy(Y_train)
            Y_test_encoded = copy.deepcopy(Y_test)
    else:
        if modeltype != 'Regression':
            mlb = My_LabelEncoder()
            if not isinstance(Y_train, str):
                Y_train_encoded= mlb.fit_transform(Y_train)
            else:
                Y_train_encoded = copy.deepcopy(Y_train)
            if not isinstance(Y_test, str):
                Y_test_encoded= mlb.transform(Y_test)
            else:
                Y_test_encoded = copy.deepcopy(Y_test)
        else:
            Y_train_encoded = copy.deepcopy(Y_train)
            Y_test_encoded = copy.deepcopy(Y_test)

    #### This is where we find out how to transform X_train and X_test  ####
    catvars = X_train.select_dtypes('object').columns.tolist() + X_train.select_dtypes('category').columns.tolist()
    numvars = X_train.select_dtypes('number').columns.tolist()
    ########    This is where we define the pipeline for cat variables ###########
    ### How do we make sure that we create one new LE_Pipe for each catvar? This is one way.
    init_str = 'make_column_transformer('
    middle_str = "".join(['(lep_one, catvars['+str(i)+']),' for i in range(len(catvars))])
    end_str = '(imp, numvars),    remainder=remainder)'
    full_str = init_str+middle_str+end_str
    ct = eval(full_str)
    pipe = make_pipeline(ct, scaler )
    ###  You will get a multidimensional numpy array ############
    dfo = pipe.fit_transform(X_train)
    if not isinstance(X_test, str):
        dfn = pipe.fit_transform(X_test)

    ### The first columns should be whatever is in the Transformer_Pipeline list of columns
    ### Hence they will be catvars. The second list will be numvars. Then only other columns that are passed through.
    ### So after the above 2 lists, you will get remainder cols unchanged: we call them leftovers.
    leftovervars = left_subtract(X_train.columns.tolist(), catvars+numvars)

    ## So if you do it correctly, you will get the list of names in proper order this way:
    ## first is catvars, then numvars and then leftovervars
    if remainder == 'drop':
        cols_names = catvars+numvars
    else:
        cols_names = catvars+numvars+leftovervars

    dfo = pd.DataFrame(dfo, columns = cols_names)
    if not isinstance(X_test, str):
        dfn = pd.DataFrame(dfn, columns = cols_names)

    copy_names = copy.deepcopy(cols_names)

    for each_col in copy_names:
        X_train_encoded[each_col] = dfo[each_col].values
        if not isinstance(X_test, str):
            X_test_encoded[each_col] = dfn[each_col].values

    return X_train_encoded, Y_train_encoded, X_test_encoded, Y_test_encoded, pipe
##################################################################################
def analyze_problem_type(train, target, verbose=0) :
    """
    ##################################################################################
    ########## Analyze if it is a Regression or Classification type problem    #######
    ##################################################################################
    """
    target = copy.deepcopy(target)
    train = copy.deepcopy(train)
    if isinstance(train, pd.Series):
        train = pd.DataFrame(train)
    cat_limit = 30 ### this determines the number of categories to name integers as classification ##
    float_limit = 15 ### this limits the number of float variable categories for it to become cat var
    if isinstance(target, str):
        target = [target]
    if len(target) == 1:
        targ = target[0]
        multilabel = False
    else:
        targ = target[0]
        multilabel = True
    ####  This is where you detect what kind of problem it is #################
    if  train[targ].dtype in ['int64', 'int32','int16']:
        if len(train[targ].unique()) <= 2:
            model_class = 'Binary_Classification'
        elif len(train[targ].unique()) > 2 and len(train[targ].unique()) <= cat_limit:
            model_class = 'Multi_Classification'
        else:
            model_class = 'Regression'
    elif  train[targ].dtype in ['float16','float32','float64']:
        if len(train[targ].unique()) <= 2:
            model_class = 'Binary_Classification'
        elif len(train[targ].unique()) > 2 and len(train[targ].unique()) <= float_limit:
            model_class = 'Multi_Classification'
        else:
            model_class = 'Regression'
    else:
        if len(train[targ].unique()) <= 2:
            model_class = 'Binary_Classification'
        else:
            model_class = 'Multi_Classification'
    ########### print this for the start of next step ###########
    if verbose >= 1:
        if multilabel:
            print('''\n########### Multi-Label %s Model Tuning and Training Started  ####''' %(model_class))
        else:
            print('''\n########### Single-Label %s Model Tuning and Training Started ####''' %(model_class))
    return model_class, multilabel
#############################################################################

reason: import numpy,import scipy,from scipy


text:
# --------------------------------------------------------
# DaSiamRPN
# Licensed under The MIT License
# Written by Qiang Wang (wangqiang2015 at ia.ac.cn)
# --------------------------------------------------------
#!/usr/bin/python

import argparse, torch.optim as optim

import cv2, torch, json
import numpy as np

import os
from os import makedirs
from os.path import realpath, dirname, join, isdir, exists

from net import SiamRPNotb
from run_SiamRPN import SiamRPN_init, SiamRPN_track
from run_SiamFPN import SiamFPN_init, SiamFPN_track
from utils import rect_2_cxy_wh, cxy_wh_2_rect
from config import *
from eval_otb import *

from fpn import *
REWRITE = False

parser = argparse.ArgumentParser(description='PyTorch SiamRPN OTB Test')
parser.add_argument('--dataset', dest='dataset', default='OTB2015', help='datasets')
parser.add_argument('-v', '--visualization', dest='visualization', action='store_true',
                    help='whether visualize result')

def train():
    return None

def track_video(model, video):
    video_path = join('test', args.dataset, 'SiamRPN_AlexNet_OTB2015')
    result_path = join(video_path, '{:s}.txt'.format(video['name']))
    if os.path.exists(result_path) and not REWRITE:
        return 1 # 跟踪过了的，就1秒咯~
    toc, regions = 0, []
    image_files, gt = video['image_files'], video['gt']
    for f, image_file in enumerate(image_files):
        im = cv2.imread(image_file)  # TODO: batch load
        tic = cv2.getTickCount()
        if f == 0:  # init
            target_pos, target_sz = rect_2_cxy_wh(gt[f])
            # state = SiamRPN_init(im, target_pos, target_sz, model)  # init tracker
            state = SiamFPN_init(im, target_pos, target_sz, model)  # init tracker
            location = cxy_wh_2_rect(state['target_pos'], state['target_sz'])
            regions.append(gt[f])
        elif f > 0:  # tracking
            # state = SiamRPN_track(state, im)  # track
            state = SiamFPN_track(state, im)  # track
            location = cxy_wh_2_rect(state['target_pos']+1, state['target_sz'])
            regions.append(location)
        toc += cv2.getTickCount() - tic

        if args.visualization and f >= 0:  # visualization
            if f == 0: cv2.destroyAllWindows()
            if len(gt[f]) == 8:
                cv2.polylines(im, [np.array(gt[f], np.int).reshape((-1, 1, 2))], True, (0, 255, 0), 3)
            else:
                cv2.rectangle(im, (gt[f, 0], gt[f, 1]), (gt[f, 0] + gt[f, 2], gt[f, 1] + gt[f, 3]), (0, 255, 0), 3)
            if len(location) == 8:
                cv2.polylines(im, [location.reshape((-1, 1, 2))], True, (0, 255, 255), 3)
            else:
                location = [int(l) for l in location]  #
                cv2.rectangle(im, (location[0], location[1]),
                              (location[0] + location[2], location[1] + location[3]), (0, 255, 255), 3)
            cv2.putText(im, str(f), (40, 40), cv2.FONT_HERSHEY_SIMPLEX, 1, (0, 255, 255), 2)

            cv2.imshow(video['name'], im)
            cv2.waitKey(1)
    toc /= cv2.getTickFrequency()

    # save result
    video_path = join('test', args.dataset, 'SiamRPN_AlexNet_OTB2015')
    if not isdir(video_path): makedirs(video_path)
    result_path = join(video_path, '{:s}.txt'.format(video['name']))
    with open(result_path, "w") as fin:
        for x in regions:
            fin.write(','.join([str(i) for i in x])+'\n')

    print('({:d}) Video: {:12s} Time: {:02.1f}s Speed: {:3.1f}fps'.format(
        v_id, video['name'], toc, f / toc))
    return f / toc

def load_dataset(dataset):
    '''[加载训练数据集，暂时以OTB100来训练]

    Arguments:
        dataset {[type]} -- [description]

    Returns:
        [type] -- [description]
    '''

    base_path = r'D:\workspace\vot\tracker_benchmark\OTB'
    json_path = r'D:\workspace\vot\tracker_benchmark\OTB\OTB2015.json'

    info = json.load(open(json_path, 'r'))
    for v in info.keys():
        path_name = info[v]['name']
        info[v]['image_files'] = [join(base_path, path_name, 'img', im_f) for im_f in info[v]['image_files']]
        info[v]['gt'] = np.array(info[v]['gt_rect'])-[1,1,0,0]  # our tracker is 0-index
        info[v]['name'] = v
    return info

def main():

    global args, v_id
    args = parser.parse_args()
    epochs = 1  # 训练的轮次
    fps_list = []
    # args.visualization = True  # 开启视觉展示

    net = SiamFPN50().cuda()

    dataset = load_dataset(args.dataset)

    optimizer = optim.SGD(net.parameters(), lr=0.001, momentum=0.9)   # 学习率为0.001
    criterion = nn.CrossEntropyLoss()   # 损失函数也可以自己定义，我们这里用的交叉熵损失函数

    for epochi in range(epochs): # 每个epoch要训练所有数据集所有的图片，每训练一个数据集便打印一下训练的效果（loss值）
        running_loss = 0.0
        for v_id, video in enumerate(dataset.keys()):
            # 此处用于测试，仅训练单个视频
            if v_id==1:
                break
            video = 'Biker'.lower()

            # 执行跟踪任务
            fps = track_video(net, dataset[video])

            # 加载跟踪后的结果
            video_path = join('test', args.dataset, 'SiamRPN_AlexNet_OTB2015')
            result_path = join(video_path, '{0}.txt'.format(video))
            bb_rect = np.array(np.loadtxt(result_path, delimiter=',').astype(np.float))

            # 加载ground_truth
            gt_rect = np.array(dataset[video]['gt_rect']).astype(np.float)

            optimizer.zero_grad()  # 梯度置零，因为反向传播过程中梯度会累加上一次循环的梯度
            loss = criterion(bb_rect, gt_rect)  # 计算损失值
            loss.backward()                     # loss反向传播
            optimizer.step()                    # 反向传播后参数更新

            print('[%d, %5d] loss: %.3f' % (epoch + 1, v_id + 1, running_loss / gt_rect.shape[0]))

            print("finish training sequence %s" % video)
            fps_list.append(fps)

            # 引入评测函数,输出精度和成功率
            # eval_auc_single(dataset='OTB2015', tracker_reg='S*', start=0, end=1e6,req_name=video)
            print('Mean Running Speed {:.1f}fps'.format(np.mean(np.array(fps_list))))
    torch.save(net.state_dict(),"siamfpn.pkl")

if __name__ == '__main__':
    main()

max_stars_repo_name: selonsy/MachineLearning

max_issues_repo_name: shenjl/MachineLearning

max_forks_repo_name: selonsy/MachineLearning


text:
import os, sys
import time
import torch
import numpy as np
sys.path.append(os.path.dirname(os.path.abspath(__file__)))
from vis_utils import get_vis_depth, get_vis_mask, get_vis_normal
import copy
import cv2
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from PIL import Image as pil
import pickle

def print_loss_pack(loss_pack, name):
    loss_depth, loss_mask_gt, loss_mask_out, loss_normal, loss_l2reg = loss_pack['depth'], loss_pack['mask_gt'], loss_pack['mask_out'], loss_pack['normal'], loss_pack['l2reg']
    if len(loss_depth.shape) == 1:
        loss_mask_gt, loss_mask_out, loss_depth, loss_normal, loss_l2reg = loss_mask_gt.mean(), loss_mask_out.mean(), loss_depth.mean(), loss_normal.mean(), loss_l2reg.mean()
    print('NAME = [{0}] -- loss_depth: {1:.4f}, loss_mask_gt: {2:.4f}, loss_mask_out: {3:.4f}, loss_normal: {4:.4f}, loss_l2reg: {5:.4f}'.format(name, loss_depth.detach().cpu().numpy(), loss_mask_gt.detach().cpu().numpy(), loss_mask_out.detach().cpu().numpy(), loss_normal.detach().cpu().numpy(), loss_l2reg.detach().cpu().numpy()))

def print_loss_pack_color(loss_pack, name):
    loss_color, loss_depth, loss_mask_gt, loss_mask_out, loss_normal, loss_l2reg, loss_l2reg_c = loss_pack['color'], loss_pack['depth'], loss_pack['mask_gt'], loss_pack['mask_out'], loss_pack['normal'], loss_pack['l2reg'], loss_pack['l2reg_c']
    print('NAME = [{0}] -- loss_color: {1:.4f}, loss_depth: {2:.4f}, loss_mask_gt: {3:.4f}, loss_mask_out: {4:.4f}, loss_normal: {5:.4f}, loss_l2reg: {6:.4f}, loss_l2re_cg: {7:.4f}'.format(name, loss_color.detach().cpu().numpy(), loss_depth.detach().cpu().numpy(), loss_mask_gt.detach().cpu().numpy(), loss_mask_out.detach().cpu().numpy(), loss_normal.detach().cpu().numpy(), loss_l2reg.detach().cpu().numpy(), loss_l2reg_c.detach().cpu().numpy()))

def demo_color_save_render_output(prefix, sdf_renderer, shape_code, color_code, camera, lighting_loc=None, profile=False):
    R, T = camera.extrinsic[:,:3], camera.extrinsic[:,3]
    R, T = torch.from_numpy(R).float().cuda(), torch.from_numpy(T).float().cuda()
    R.requires_grad, T.requires_grad = False, False

    if lighting_loc is not None:
        lighting_locations = torch.from_numpy(lighting_loc).float().unsqueeze(0).cuda()
    else:
        lighting_locations = None
    render_output = sdf_renderer.render(color_code, shape_code, R, T, profile=profile, no_grad=True, lighting_locations=lighting_locations)
    depth_rendered, normal_rendered, color_rgb, valid_mask_rendered, min_sdf_sample = render_output

    data = {}
    data['depth'] = depth_rendered.detach().cpu().numpy()
    data['normal'] = normal_rendered.detach().cpu().numpy()
    data['mask'] = valid_mask_rendered.detach().cpu().numpy()
    data['color'] = color_rgb.detach().cpu().numpy()
    data['min_sdf_sample'] = min_sdf_sample.detach().cpu().numpy()
    data['latent_tensor'] = shape_code.detach().cpu().numpy()
    data['K'] = sdf_renderer.get_intrinsic()
    data['RT'] = torch.cat([R, T[:,None]], 1).detach().cpu().numpy()
    fname = prefix + '_info.pkl'
    with open(fname, 'wb') as f:
        pickle.dump(data, f)

    img_hw = sdf_renderer.get_img_hw()
    visualizer = Visualizer(img_hw)
    print('Writing to prefix: {}'.format(prefix))
    visualizer.visualize_depth(prefix + '_depth.png', depth_rendered.detach().cpu().numpy(), valid_mask_rendered.detach().cpu().numpy())
    visualizer.visualize_normal(prefix + '_normal.png', normal_rendered.detach().cpu().numpy(), valid_mask_rendered.detach().cpu().numpy(), bgr2rgb=True)
    visualizer.visualize_mask(prefix + '_silhouette.png', valid_mask_rendered.detach().cpu().numpy())
    cv2.imwrite(prefix + '_rendered_rgb.png', color_rgb.detach().cpu().numpy() * 255)

class Visualizer(object):
    def __init__(self, img_hw, dmin=0.0, dmax=10.0):
        self.img_h, self.img_w = img_hw[0], img_hw[1]
        self.data = {}
        self.dmin, self.dmax = dmin, dmax

        self.loss_counter = 0
        self.loss_curve = {}
        self.loss_list = []
        self.chamfer_list = []

    def get_data(self, data_name):
        if data_name in self.data.keys():
            return self.data[data_name]
        else:
            raise ValueError('Key {0} does not exist.'.format(data_name))

    def set_data(self, data):
        self.data = data

    def reset_data(self):
        self.data = {}
        keys = ['mask_gt', 'mask_output', 'loss_mask_gt', 'loss_mask_out',
                'depth_gt', 'depth_output', 'loss_depth',
                'normal_gt', 'normal_output', 'loss_normal']
        for key in keys:
            self.data[key] = np.zeros((64, 64))

    def reset_loss_curve(self):
        self.loss_counter = 0
        self.loss_curve = {}

    def reset_all(self):
        self.reset_data()
        self.reset_loss_curve()

    def add_loss_from_pack(self, loss_pack):
        '''
        potential properties:
        ['mask_gt', 'mask_out', 'depth' 'normal', 'l2reg']
        '''
        loss_name_list = list(loss_pack.keys())
        if self.loss_curve == {}:
            for loss_name in loss_name_list:
                self.loss_curve[loss_name] = []
        for loss_name in loss_name_list:
            loss_value = loss_pack[loss_name].detach().cpu().numpy()
            self.loss_curve[loss_name].append(loss_value)
        self.loss_counter = self.loss_counter + 1

    def add_loss(self, loss):
        self.loss_list.append(loss.detach().cpu().numpy())

    def add_chamfer(self, chamfer):
        self.chamfer_list.append(chamfer)

    def add_data(self, data_name, data_src, data_mask=None):
        '''
        potential properties:
        mask: 	['mask_gt', 'mask_output', 'loss_mask_gt', 'loss_mask_out']
        depth: 	['depth_gt', 'depth_output', 'loss_depth']
        normal:	['normal_gt', 'normal_output', 'loss_normal']
        '''
        if data_mask is None:
            self.data[data_name] = data_src
        else:
            data_map = np.zeros(data_mask.shape)
            data_map[data_mask != 0] = data_src
            self.data[data_name] = data_map

    def save_depth(self, fname, depth_vis, cmap='magma', direct=False):
        if direct:
            cv2.imwrite(fname, depth_vis)
            return 0
        vmin, vmax = 0, 255
        normalizer = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
        mapper = cm.ScalarMappable(norm=normalizer, cmap=cmap)
        colormapped_im = (mapper.to_rgba(depth_vis)[:,:,:3] * 255).astype(np.uint8)
        im = pil.fromarray(colormapped_im)
        im.save(fname)

    def save_mask(self, fname, mask_vis, bgr2rgb=False):
        if bgr2rgb:
            mask_vis = cv2.cvtColor(mask_vis, cv2.COLOR_BGR2RGB)
        cv2.imwrite(fname, mask_vis)

    def save_normal(self, fname, normal_vis, bgr2rgb=False):
        if bgr2rgb:
            normal_vis = cv2.cvtColor(normal_vis, cv2.COLOR_BGR2RGB)
        cv2.imwrite(fname, normal_vis)

    def save_error(self, fname, error_vis, bgr2rgb=False):
        self.save_depth(fname, error_vis, cmap='jet')

    def visualize_depth(self, fname, depth, mask=None):
        # depth_vis = get_vis_depth(depth, mask=mask, dmin=self.dmin, dmax=self.dmax)
        depth_vis = get_vis_depth(depth, mask=mask)
        # self.save_depth(fname, depth_vis)
        cv2.imwrite(fname, depth_vis)

    def visualize_normal(self, fname, normal, mask=None, bgr2rgb=False):
        normal_vis = get_vis_normal(normal, mask=mask)
        if bgr2rgb:
            normal_vis = cv2.cvtColor(normal_vis, cv2.COLOR_BGR2RGB)
        cv2.imwrite(fname, normal_vis)

    def visualize_mask(self, fname, mask, bgr2rgb=False):
        mask_vis = get_vis_mask(mask)
        if bgr2rgb:
            mask_vis = cv2.cvtColor(mask_vis, cv2.COLOR_BGR2RGB)
        cv2.imwrite(fname, mask_vis)

    def imshow(self, ax, img, title=None):
        ax.imshow(img)
        ax.axis('off')
        if title is not None:
            ax.set_title(title)

    def imshow_bgr2rgb(self, ax, img, title=None):
        if len(img.shape) == 3:
            img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        ax.imshow(img)
        ax.axis('off')
        if title is not None:
            ax.set_title(title)

    def show_loss_curve(self, fname):
        pass

    def show_all_data_3x4(self, fname):
        fig, axs = plt.subplots(3, 4, figsize=(30,30))

        # first row, groundtruth
        depth_gt_vis = get_vis_depth(self.data['depth_gt'], mask=self.data['mask_gt'], dmin=self.dmin, dmax=self.dmax)
        self.imshow_bgr2rgb(axs[0, 0], 255 - depth_gt_vis, title='depth gt')
        normal_gt_vis = get_vis_normal(self.data['normal_gt'], mask=self.data['mask_gt'])
        self.imshow(axs[0, 1], normal_gt_vis, title='normal gt')
        mask_gt_vis = get_vis_mask(self.data['mask_gt'])
        self.imshow_bgr2rgb(axs[0, 2], 255 - mask_gt_vis, title='mask gt')
        axs[0, 3].axis('off')

        # second row, output
        depth_output_vis = get_vis_depth(self.data['depth_output'], mask=self.data['mask_output'], dmin=self.dmin, dmax=self.dmax)
        self.imshow_bgr2rgb(axs[1, 0], 255 - depth_output_vis, title='depth output')
        normal_output_vis = get_vis_normal(self.data['normal_output'], mask=self.data['mask_output'])
        self.imshow(axs[1, 1], normal_output_vis, title='normal output')
        mask_output_vis = get_vis_mask(self.data['mask_output'])
        self.imshow_bgr2rgb(axs[1, 2], 255 - mask_output_vis, title='mask output')
        axs[1, 3].axis('off')

        # third row, loss
        valid_mask = np.logical_and(self.data['mask_gt'], self.data['mask_output'])
        loss_depth_vis = get_vis_depth(np.abs(self.data['loss_depth']), valid_mask, dmin=0.0, dmax=0.5)
        self.imshow_bgr2rgb(axs[2, 0], 255 - loss_depth_vis, title='depth loss')
        loss_normal_vis = get_vis_depth(self.data['loss_normal'], valid_mask, dmin=-1.0, dmax=0.0)
        self.imshow_bgr2rgb(axs[2, 1], 255 - loss_normal_vis, title='normal loss')
        loss_mask_gt_vis = get_vis_mask(np.abs(self.data['loss_mask_gt']) > 0)
        self.imshow_bgr2rgb(axs[2, 2], 255 - loss_mask_gt_vis, title='gt \ output')
        loss_mask_out_vis = get_vis_mask(np.abs(self.data['loss_mask_out']) > 0)
        self.imshow_bgr2rgb(axs[2, 3], 255 - loss_mask_out_vis, title='output \ gt')

        # savefig
        fig.savefig(fname)
        plt.close('all')

    def save_all_data(self, prefix):
        # groundtruth
        depth_gt_vis = get_vis_depth(self.data['depth_gt'], mask=self.data['mask_gt'], dmin=self.dmin, dmax=self.dmax)
        self.save_depth(prefix + '_depth_gt.png', depth_gt_vis, cmap='magma', direct=True)
        normal_gt_vis = get_vis_normal(self.data['normal_gt'], mask=self.data['mask_gt'])
        self.save_normal(prefix + '_normal_gt.png', normal_gt_vis, bgr2rgb=True)
        mask_gt_vis = get_vis_mask(self.data['mask_gt'])
        self.save_mask(prefix + '_mask_gt.png', mask_gt_vis)

        # output
        depth_output_vis = get_vis_depth(self.data['depth_output'], mask=self.data['mask_output'], dmin=self.dmin, dmax=self.dmax)
        self.save_depth(prefix + '_depth_output.png', depth_output_vis, cmap='magma', direct=True)
        normal_output_vis = get_vis_normal(self.data['normal_output'], mask=self.data['mask_output'])
        self.save_normal(prefix + '_normal_output.png', normal_output_vis, bgr2rgb=True)
        mask_output_vis = get_vis_mask(self.data['mask_output'])
        self.save_mask(prefix + '_mask_output.png', mask_output_vis)

        # third row, loss
        valid_mask = np.logical_and(self.data['mask_gt'], self.data['mask_output'])
        loss_depth_vis = get_vis_depth(np.abs(self.data['loss_depth']), valid_mask, dmin=0.0, dmax=0.5, bg_color=0)
        self.save_error(prefix + '_depth_loss.png', loss_depth_vis, bgr2rgb=True)
        loss_normal_vis = get_vis_depth(self.data['loss_normal'], valid_mask, dmin=-1.0, dmax=0.0, bg_color=0)
        self.save_error(prefix + '_normal_loss.png', loss_normal_vis, bgr2rgb=True)

        loss_mask_gt_vis = get_vis_depth(np.abs(self.data['loss_mask_gt']), bg_color=0)
        self.save_error(prefix + '_mask_gt_loss.png', loss_mask_gt_vis, bgr2rgb=True)
        loss_mask_out_vis = get_vis_depth(np.abs(self.data['loss_mask_out']), bg_color=0)
        self.save_error(prefix + '_mask_out_loss.png', loss_mask_out_vis, bgr2rgb=True)
        self.save_error(prefix + '_mask_loss.png', loss_mask_gt_vis + loss_mask_out_vis, bgr2rgb=True)

    def dump_all_data(self, fname):
        with open(fname, 'wb') as f:
            pickle.dump({'data': self.data, 'loss_curve': self.loss_curve, 'loss_list': self.loss_list, 'chamfer_list': self.chamfer_list}, f)

    def show_all_data(self, fname):
        self.show_all_data_3x4(fname)
        # self.save_all_data(fname[:-4])

    def show_all_data_color(self, fname):
        fig, axs = plt.subplots(3, 4, figsize=(30,30))

        # first row, groundtruth
        depth_gt_vis = get_vis_depth(self.data['depth_gt'], mask=self.data['mask_gt'], dmin=self.dmin, dmax=self.dmax)
        self.imshow_bgr2rgb(axs[0, 0], depth_gt_vis, title='depth gt')
        normal_gt_vis = get_vis_normal(self.data['normal_gt'])
        self.imshow_bgr2rgb(axs[0, 1], normal_gt_vis, title='normal gt')
        mask_gt_vis = get_vis_mask(self.data['mask_gt'])
        self.imshow_bgr2rgb(axs[0, 2], mask_gt_vis, title='mask gt')
        self.imshow_bgr2rgb(axs[0, 3], self.data['color_gt'], title='rgb gt')

        # second row, output
        depth_output_vis = get_vis_depth(self.data['depth_output'], mask=self.data['mask_output'], dmin=self.dmin, dmax=self.dmax)
        self.imshow_bgr2rgb(axs[1, 0], depth_output_vis, title='depth output')
        normal_output_vis = get_vis_normal(self.data['normal_output'])
        self.imshow_bgr2rgb(axs[1, 1], normal_output_vis, title='normal output')
        mask_output_vis = get_vis_mask(self.data['mask_output'])
        self.imshow_bgr2rgb(axs[1, 2], mask_output_vis, title='mask output')
        self.imshow_bgr2rgb(axs[1, 3], self.data['color_output'], title='rgb output')

        # third row, loss
        valid_mask = np.logical_and(self.data['mask_gt'], self.data['mask_output'])
        loss_depth_vis = get_vis_depth(np.abs(self.data['loss_depth']), valid_mask, dmin=0.0, dmax=0.5)
        self.imshow_bgr2rgb(axs[2, 0], loss_depth_vis, title='depth loss')
        loss_normal_vis = get_vis_depth(self.data['loss_normal'], valid_mask, dmin=-1.0, dmax=0.0)
        self.imshow_bgr2rgb(axs[2, 1], loss_normal_vis, title='normal loss')
        loss_mask_gt_vis = get_vis_mask(np.abs(self.data['loss_mask_gt']) > 0)
        loss_mask_out_vis = get_vis_mask(np.abs(self.data['loss_mask_out']) > 0)
        loss_mask_gt_vis += loss_mask_out_vis
        self.imshow_bgr2rgb(axs[2, 2], loss_mask_gt_vis, title='mask loss')
        self.imshow_bgr2rgb(axs[2, 3], self.data['loss_color'], title='rgb loss')

        # savefig
        fig.savefig(fname)
        plt.close('all')

    def return_output_data_color(self):
        return self.data['color_output'], self.data['depth_output'], self.data['normal_output'], self.data['mask_output']

    def show_all_data_color_multi(self, fname, num_img=4):
        fig, axs = plt.subplots(3, 2*num_img, figsize=(8*2*num_img,25))

        for i in range(num_img):
            # first row, ground truth
            self.imshow_bgr2rgb(axs[0, 2*i], self.data['color_gt-{}'.format(i)], title='rgb gt {}'.format(i))
            mask_gt_vis = get_vis_mask(self.data['mask_gt-{}'.format(i)])
            self.imshow_bgr2rgb(axs[0, 2*i+1], mask_gt_vis, title='mask gt {}'.format(i))

            # second row, output
            self.imshow_bgr2rgb(axs[1, 2*i], self.data['color_output-{}'.format(i)], title='rgb output {}'.format(i))
            mask_output_vis = get_vis_mask(self.data['mask_output-{}'.format(i)])
            self.imshow_bgr2rgb(axs[1, 2*i+1], mask_output_vis, title='mask output {}'.format(i))

            # third row, loss
            self.imshow_bgr2rgb(axs[2, 2*i], self.data['loss_color-{}'.format(i)], title='rgb loss {}'.format(i))
            loss_mask_gt_vis = get_vis_mask(np.abs(self.data['loss_mask_gt-{}'.format(i)]) > 0)
            loss_mask_out_vis = get_vis_mask(np.abs(self.data['loss_mask_out-{}'.format(i)]) > 0)
            loss_mask_gt_vis += loss_mask_out_vis
            self.imshow_bgr2rgb(axs[2, 2*i+1], loss_mask_gt_vis, title='mask loss {}'.format(i))

        # savefig
        plt.subplots_adjust(top=0.95, right=0.99, left=0.01, bottom=0.01, wspace=0.05, hspace=0.1)
        fig.savefig(fname)
        plt.close('all')

    def show_all_data_color_warp(self, fname):
        fig, axs = plt.subplots(1, 5, figsize=(15, 3.4))
        self.imshow_bgr2rgb(axs[0], self.data['color_gt-1'], title='view 1')
        self.imshow_bgr2rgb(axs[1], self.data['color_gt-2'], title='view 2')
        self.imshow_bgr2rgb(axs[2], self.data['color_valid-1'], title='valid region in view 1')
        self.imshow_bgr2rgb(axs[3], self.data['color_valid-2'], title='warped color from view 2')
        self.imshow_bgr2rgb(axs[4], self.data['color_valid_loss'], title='color loss')

        # savefig
        plt.subplots_adjust(top=0.99, right=0.99, left=0.01, bottom=0.00, wspace=0.05, hspace=0)
        fig.savefig(fname)
        plt.close('all')

max_stars_repo_name: hyunynim/DIST-Renderer

max_issues_repo_name: hyunynim/DIST-Renderer

max_forks_repo_name: hyunynim/DIST-Renderer


text:
# Script to statistically asess the industrial potential of heterologous reactions
# Author: Szabolcs Cselgo Kovacs
# Contact: kovszasz@gmail.com

import sys,cobra, cobra.io, cobra.test, csv,pandas,codecs, email, smtplib,sys,random,pickle
#from cameo import *
#from cameo import load_model
from cobra import Reaction, Metabolite, Model
import cobra, cobra.io, cobra.test, csv,pandas,codecs,pickle
from cobra.util.solver import linear_reaction_coefficients
from random import randint
#import cameo
import numpy as np
#model=load_model("heterologous.xml")
model="../heterologous.xml"
model=cobra.io.read_sbml_model(model)
model.solver='gurobi'

def white_space_remover(string):#while reading the '.txt' files, the lines are separated wit the '\n' character. This function removes '\n' from the end of the strings.
    return string.replace('\n','')

def binary_saver(obj,name):#We store the results and other important variables binary. It is easer and faster to read and save.
    f=''
    outfile=open(f+name,'wb')
    pickle.dump(obj,outfile)
    outfile.close()

def binary_loader(name):#We store the results and other important variables binary. It is easer and faster to read and save.
    f=''
    infile=open(f+name,'rb')
    return pickle.load(infile)

heterologous_reactions={}#the key is the ID of the heterologous reaction, the item is information about its reversibility.
for i in model.reactions:
    if i.id.find('Meta')>-1:
        heterologous_reactions[i.id]=i.reversibility
        i.knock_out() #all heterologous reaction will be knocked out, because the randomly generated reaction sets are loaded by the path_random input

path_random='RandomHeterologousReactionSet' #input("Reaction set:") #a path for the randomly generated heterologous reaction sets
random_reaction_set=binary_loader(str(path_random)) #loading the randomly generated heterologous reaction sets
class ProductionPrediction():
    def __init__(self,model,target,C_source,random_set,start_seed,end_seed,heterologous_reactions):
        """
        model: heterologous metabolic network
        target: the ID of the target compound
        random_set: the randomly generated heterologous reaction set
        start_seed: the random seeds are between 8011-9011 (1000 randomly generated reaction set)
        end_seed: the random seeds are between 8011-9011 (1000 randomly generated reaction set)
        heterologous_reactions: dict of the heterologous reactions, containing information about their reversibility
        """
        self.model=model
        self.model.reactions.get_by_id("Ec_biomass_iJO1366_core_53p95M").lower_bound=0.09915628719962316 #10% of the maximal biomass production value of native metabolic network. This value based on the compound_production.py results
        self.model.reactions.get_by_id(C_source).lower_bound=-10
        self.model.reactions.get_by_id(C_source).upper_bound=1000
        self.heterologous_reactions=heterologous_reactions
        self.target=target
        self.C_source=C_source
        self.seed_list=[i for i in range(int(start_seed),int(end_seed))] #generating the random seed range
        self.reaction_set=random_set #a dictionary, where the key is a random seed number, and the item is a list containing the randomly generated heterologous reactions
        self.result=[]
        self.create_target_demand_reaction(self.target)
        for i in self.seed_list: #iterating the among the random seed values
            self.result.append(self.run_prediction(self.target,self.C_source,self.reaction_set[i]))
        self.model.remove_reactions('DEMAND_target')

    def c_atom_number(self,a):#For the yield calculation it is required to get information about the carbon atom numbe of the given compund. This function calculates it from the formula of the compound.
        C={}
        for i in ['Cl','H','O','S','P','C','N']:
            if a.find(i)>-1:
                a=a.replace(i,','+i+';')
        l1=a.split(',')
        for i in l1:
            p=i.split(';')
            if len(p)==2:
                if p[1]=='':
                    C[p[0]]=1
                else:
                    C[p[0]]=int(p[1])
        return C['C']

    def create_target_demand_reaction(self,target):
        target_metabol=self.model.metabolites.get_by_id(target)
        target_demand=Reaction(id="DEMAND_target")
        target_demand.add_metabolites({target_metabol:-1.0})
        self.model.add_reaction(target_demand)
        self.model.reactions.get_by_id("DEMAND_target").lower_bound=0
        self.model.reactions.get_by_id("DEMAND_target").upper_bound=1000
        self.model.objective={self.model.reactions.DEMAND_target:1}

    def run_prediction(self,target,C_source,heterologous_reaction_list):
        """
        This function simulates the production of the target compound, while only the randomly generated heterologous reaction set is turned on
        """
        for i in self.heterologous_reactions.keys(): # if any undesired heterologous reaction might be turned on (undesired: not part of the randomly generated heterologous reaction list)
            self.model.reactions.get_by_id(i).knock_out()
        for i in heterologous_reaction_list:
            if self.heterologous_reactions[i]:
                self.model.reactions.get_by_id(i).lower_bound=-1000
                self.model.reactions.get_by_id(i).upper_bound=1000
            else:
                self.model.reactions.get_by_id(i).lower_bound=0
                self.model.reactions.get_by_id(i).upper_bound=1000
        try:
            p=cobra.flux_analysis.pfba(self.model)

            print(p.fluxes['DEMAND_target'])
            compound_yield=abs((p.fluxes["DEMAND_target"]/p.fluxes[C_source])*(self.c_atom_number(self.model.metabolites.get_by_id(target).formula)/6)) # 6 represents the carbon atom number of glucose (default carbon soruce)
        except:
            compound_yield=0
        return compound_yield
start_seed=sys.argv[1]
end_seed=sys.argv[2]
name=sys.argv[3]

targets=[]
with open("value_added_compound_list_for_heterologous_model.txt") as target_compounds:
    for i in target_compounds:
        targets.append(white_space_remover(i))

result={}
#targets=['ptrc_c']
carbon_source="EX_glc_lp_e_rp_"
for target in targets:
    print(target)
    result[target]=ProductionPrediction(model,target,carbon_source,random_reaction_set,start_seed,end_seed,heterologous_reactions).result

binary_saver(result,'results/'+str(name))


text:
import numpy as np
#import RPi.GPIO as gpio
import needle_detection.parameters as p
from skimage import transform
import matplotlib.pyplot as plt
import cv2 as cv

#def get_angle_from_pi():
#    gpio.setmode(gpio.BOARD)
#    gpio.setup(p.channels, gpio.IN, pull_up_down=gpio.PUD_UP)
#    gpio.add_event_detect(p.channels[0], gpio.FALLING,  bouncetime=200)
#    print("Set up input channels")
#    while True:

     # Angle change button pressed ?
#     if gpio.event_detected(p.channels[0]):
#         print("Button pressed!")
#         # Check each input
#         for c in range(2, len(p.channels)):
#             # If input is low the expected angle is updated
#             if gpio.input(p.channels[c]) == gpio.LOW:
#                 expected_angle = p.angles[c - 1]
#                 insertion_depth = p.insertion_depths[c - 1]
#                 print("Angle set to ", expected_angle)
#                 new_angle = 1
#                 break
#    return

#def read_cam(channel=0):
#    cap = cv2.VideoCapture(channel)
#    cap.open(channel)
#    cap.set(cv2.CAP_PROP_FRAME_WIDTH, p.image_size_x)
#    cap.set(cv2.CAP_PROP_FRAME_HEIGHT, p.image_size_y)
#    print("Capture device opened!")#

                # Get image and crop to ROI
#            ret, frame = cap.read()
#            if ret == False:
#                print("Video grabber disconnected!")
#                exit
#    return

#uhrzeigersinn, startet oben links
def mask_roi(img, vertices):
    """
    Get roi using a mask

    Parameters
    ----------
    img : numpy.ndarray
        Frame of us image
    vertices: numpy.ndarray
        Vertices characterize the polynomial, defined in parameters.py
    Returns
    -------
    numpy.ndarray
        masked input image
    """
    #blank mask:
    mask = np.zeros_like(img)

    # fill the mask
    #cv.fillPoly(mask, vertices, 255)
    cv.fillConvexPoly(mask, np.array(vertices, 'int32'), 255)
    # now only show the area that is the mask
    masked = cv.bitwise_and(img, mask)

    return masked

def read_frame(path, alpha=0, beta=255):
    """
    Reads in and normalizes single frame in range(0,255)

    Parameters
    ----------
    path : str
        The file location of the frame

    Returns
    -------
    numpy.ndarray
        normalized gray-scale image
    """

    frame = cv.imread(path, cv.IMREAD_GRAYSCALE)
    norm = cv.normalize(frame, None, alpha, beta, norm_type=cv.NORM_MINMAX, dtype=cv.CV_32F)

    return norm

#def roi(frame, vertices):
    #blank mask:
 #   mask = np.zeros_like(img)
    # fill the mask
 #   cv.fillPoly(mask, vertices, 255)
    # now only show the area that is the mask
#    masked = cv.bitwise_and(img, mask)
#    return masked

def get_ROI(frame, angle):
    """
    Defines ROI in given frame

    Parameters
    ----------
    frame : numpy.ndarray
    angle: int
         Angel in degrees of the needle holder measuered with respect to 'vertical' transducer axis

    Returns
    -------
    frame_roi : numpy.ndarray
        ROI of given frame
    rescale_factor: float
        xxx
    """

    width = p.roi_x_2 - p.roi_x_1
    height = p.roi_y_2 - p.roi_y_1

    dist_per_pixel = p.depth/height     # Distance in cm per Pixel
    angle = np.pi/2-np.deg2rad(angle)   # Recalc expected angle so it is measured with respect to the x-axis

    # Calculate new image size
    if width <= height:
        new_height = p.processing_size
        rescale_factor = new_height/height
        new_width = round(rescale_factor*width)
    else:
        new_width = p.processing_size
        rescale_factor = new_width/width
        new_height = round(rescale_factor*height)
    width = new_width
    height = new_height

    wdt_hgt_ref = np.sqrt(width**2+height**2)
    frame_raw = frame[p.roi_y_1:p.roi_y_2, p.roi_x_1:p.roi_x_2]
    frame_roi = transform.resize(frame_raw, (height, width))

    return frame_roi, rescale_factor