text:
# coding: utf-8

# In[ ]:

from keras.layers import Input, Dense, Flatten, Dropout
from keras.layers.core import Reshape
from keras.models import Model
from keras.callbacks import ModelCheckpoint
from keras.layers.convolutional import MaxPooling2D,UpSampling2D,Conv2DTranspose
from keras.layers.convolutional import Convolution2D as Conv2D
from keras.layers.normalization import BatchNormalization
from keras.optimizers import Adadelta, RMSprop, adam, SGD
from keras import backend as K
from keras.models import model_from_json
from keras.utils import np_utils
import os
import os.path
import numpy as np
from PIL import Image
from numpy import *
from sklearn.utils import shuffle
from sklearn.cross_validation import train_test_split
import scipy.misc
import cv2
import matplotlib.pyplot as plt
import matplotlib
#from sklearn.metrics import classification_report,confusion_matrix

# In[ ]:

input_img = Input(shape=(128,128,1))
num_classes = 123
img_rows, img_cols = 128, 128

# In[ ]:

enco = Conv2D(16, (3, 3), activation='relu', padding='same')(input_img)
enco = BatchNormalization()(enco)
enco = Conv2D(16, (3, 3), activation='relu', padding='same')(enco)
enco = BatchNormalization()(enco)
enco = MaxPooling2D(pool_size=(2, 2))(enco)

enco = Conv2D(32, (3, 3), activation='relu', padding='same')(enco)
enco = BatchNormalization()(enco)
enco = Conv2D(32, (3, 3), activation='relu', padding='same')(enco)
enco = BatchNormalization()(enco)
enco = MaxPooling2D(pool_size=(2, 2))(enco)

enco = Conv2D(64, (3, 3), activation='relu', padding='same')(enco)
enco = BatchNormalization()(enco)
enco = Conv2D(64, (3, 3), activation='relu', padding='same')(enco)
enco = BatchNormalization()(enco)
enco = MaxPooling2D(pool_size=(2, 2))(enco)

enco = Conv2D(128, (3, 3), activation='relu', padding='same')(enco)
enco = BatchNormalization()(enco)
enco = Conv2D(128, (3, 3), activation='relu', padding='same')(enco)
enco = BatchNormalization()(enco)

#enco.load_weights("Only_Encoder.h5")
encoder = Model(input_img, enco)
encoder.load_weights("Only_Encoder.h5")

classify = Flatten()(enco)
classify = Dense(64, activation='relu')(classify)
classify = Dense(32, activation='relu')(classify)
classify = Dense(num_classes, activation='softmax')(classify)

#network = Model(enco, classify)
network = Model(input_img, classify)
rms=RMSprop(lr=0.001, rho=0.9, epsilon=1e-08, decay=0.001)
#network.compile(loss='mean_squared_error', optimizer=rms)
network.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])
#network.summary()
for layers in encoder.layers:
	layers.trainable=False

#network.summary()
#exit()

# In[ ]:

path="Data"
basic_mat=[]
index=[]
epochs=100
batch_size1=64

# In[ ]:

for i in range(1,124):
    path_major=path+'/'+str(i)
    for j in range(1,101):
        img=array(Image.open(path_major+"/"+str(j)+"_.jpg"))
        #print shape(img)
        img = cv2.cvtColor( img, cv2.COLOR_RGB2GRAY )
        img=img.reshape(128,128,1)
        basic_mat.append(img)
        index.append(i-1)

network.summary()

# In[ ]:

data,Label = shuffle(basic_mat,index, random_state=2)

# In[ ]:

X_train, X_test, y_train, y_test = train_test_split(data, Label, test_size=0.2, random_state=2)
X_train = array(X_train)
y_train = array(y_train)
X_test = array(X_test)
y_test = array(y_test)
print(X_train.shape)
print(y_train.shape)
print(X_test.shape)
print(y_test.shape)

y_train = np_utils.to_categorical(y_train)
y_test = np_utils.to_categorical(y_test)

# In[ ]:

x_train = X_train.astype('float32') / 255.
x_test = X_test.astype('float32') / 255.
network.fit(x_train, y_train, validation_data=(x_test, y_test),batch_size=batch_size1, nb_epoch=epochs, verbose=2)
scores = network.evaluate(x_test, y_test, verbose=0)
print ("%s: %.2f%%" % (network.metrics_names[1], scores[1]*100))
network.save_weights("only_classify.h5")
# final_network.fit(x_train, y_train, validation_data=(x_test, y_test),batch_size=batch_size1, nb_epoch=epochs, verbose=2)
# scores_final = final_network.evaluate(x_test, y_test, verbose=0)
# print ("%s: %.2f%%" % (final_network.metrics_names[1], scores_final[1]*100))
# final_network.save_weights("complete_model.h5")

max_stars_repo_name: shreyasbapat/Fabric-Detection

max_issues_repo_name: shreyasbapat/Fabric-Detection

max_forks_repo_name: shreyasbapat/Fabric-Detection

reason: import numpy,from numpy,import scipy


text:
import xlwings as xw
import time
import numpy as np
import pandas as pd

@xw.func(async_mode='threading')
@xw.arg("n", numbers=int)
@xw.arg("m", numbers=int)
@xw.ret(expand='table')
def threading(n, m):
    print('1 - CALLED THREADING FUNCTION')
    time.sleep(1)
    return [
        ["%s x %s : %s, %s" % (n, m, i, j) for j in range(m)]
        for i in range(n)
    ]

@xw.arg("n", numbers=int)
@xw.arg("m", numbers=int)
@xw.ret(expand='table')
def nothreading(n, m):
    print('2 - CALLED NOTHREADING FUNCTION')
    return [
        ["%s x %s : %s, %s" % (n, m, i, j) for j in range(m)]
        for i in range(n)
    ]

@xw.func
def simple(x):
    print('3 - CALLED SIMPLE FUNCTION')
    return x

@xw.func(async_mode='threading')
def simple_threading(x):
    print('4 - CALLED SIMPLE THREADING FUNCTION')
    time.sleep(1)
    return x

@xw.func(async_mode='threading')
@xw.arg("i", numbers=int)
@xw.arg("j", numbers=int)
@xw.ret(expand='table')
def numpy_async(i, j):
    print('5 - CALLED NUMPY ASYNC FUNCTION')
    time.sleep(1)
    return np.arange(i * j).reshape((i, j))

if __name__ == '__main__':
    xw.serve()

ext: py

max_stars_repo_name: kushal-kumaran/xlwings

max_issues_repo_name: kushal-kumaran/xlwings

max_forks_repo_name: kushal-kumaran/xlwings

reason: import numpy


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

FONTSIZE = 18

def plot_dynamics(history_pi_0, history_pi_1, pi_alpha_gradient_history, pi_beta_gradient_history, title=''):
    cmap = plt.get_cmap('viridis')
    colors = range(len(history_pi_1))
    fig = plt.figure(figsize=(6, 5))
    ax = fig.add_subplot(111)

    scatter = ax.scatter(history_pi_0, history_pi_1, c=colors, s=1)
    ax.scatter(0.5, 0.5, c='r', s=10., marker='*')
    colorbar = fig.colorbar(scatter, ax=ax)
    colorbar.set_label('Iterations', rotation=270, fontsize=FONTSIZE)

    skip = slice(0, len(history_pi_0), 5)
    ax.quiver(history_pi_0[skip],
              history_pi_1[skip],
              pi_alpha_gradient_history[skip],
              pi_beta_gradient_history[skip],
              units='xy', scale=10., zorder=3, color='blue',
              width=0.007, headwidth=3., headlength=4.)

    ax.set_ylabel("Policy of Player 2", fontsize=FONTSIZE)
    ax.set_xlabel("Policy of Player 1", fontsize=FONTSIZE)
    ax.set_ylim(0, 1)
    ax.set_xlim(0, 1)
    ax.set_xticklabels(map(str, np.arange(0, 11, 2) / 10), fontsize=FONTSIZE-6)
    ax.set_yticklabels(map(str, np.arange(0, 11, 2) / 10), fontsize=FONTSIZE-6)
    ax.set_title(title, fontsize=FONTSIZE+8)
    plt.tight_layout()
    plt.savefig('{}.pdf'.format(title))
    plt.show()

def plot_Q(data, row_labels, col_labels, title, ax=None,
          cbar_kw={}, cbarlabel="", **kwargs):
    fig = plt.figure(figsize=(6.5, 5))
    ax = fig.add_subplot(111)

    im, cbar = heatmap(data, row_labels, col_labels, ax=ax,
                       cmap="YlGn", cbarlabel="Q Values")
    texts = annotate_heatmap(im, valfmt="{x:.1f} t")
    # print(data)
    ax.set_title(title, pad=50)
    ax.xaxis.set_label_position("top")
    ax.set_xlabel('Alpha')
    ax.set_ylabel('Beta')
    fig.tight_layout()
    plt.savefig('{}.pdf'.format(title), bbox_inches='tight', format='pdf')
    plt.show()

def heatmap(data, row_labels, col_labels, ax=None,
            cbar_kw={}, cbarlabel="", **kwargs):
    """
    Create a heatmap from a numpy array and two lists of labels.

    Arguments:
        data       : A 2D numpy array of shape (N,M)
        row_labels : A list or array of length N with the labels
                     for the rows
        col_labels : A list or array of length M with the labels
                     for the columns
    Optional arguments:
        ax         : A matplotlib.axes.Axes instance to which the heatmap
                     is plotted. If not provided, use current axes or
                     create a new one.
        cbar_kw    : A dictionary with arguments to
                     :meth:`matplotlib.Figure.colorbar`.
        cbarlabel  : The label for the colorbar
    All other arguments are directly passed on to the imshow call.
    """

    if not ax:
        ax = plt.gca()

    # Plot the heatmap
    im = ax.imshow(data, **kwargs)

    # Create colorbar
    cbar = ax.figure.colorbar(im, ax=ax, **cbar_kw)
    cbar.ax.set_ylabel(cbarlabel, rotation=-90, va="bottom")

    # We want to show all ticks...
    ax.set_xticks(np.arange(data.shape[1]))
    ax.set_yticks(np.arange(data.shape[0]))
    # ... and label them with the respective list entries.
    ax.set_xticklabels(col_labels)
    ax.set_yticklabels(row_labels)

    # Let the horizontal axes labeling appear on top.
    ax.tick_params(top=True, bottom=False,
                   labeltop=True, labelbottom=False)

    # Rotate the tick labels and set their alignment.
    plt.setp(ax.get_xticklabels(), rotation=-30, ha="right",
             rotation_mode="anchor")

    # Turn spines off and create white grid.
    for edge, spine in ax.spines.items():
        spine.set_visible(False)

    ax.set_xticks(np.arange(data.shape[1]+1)-.5, minor=True)
    ax.set_yticks(np.arange(data.shape[0]+1)-.5, minor=True)
    ax.grid(which="minor", color="w", linestyle='-', linewidth=3)
    ax.tick_params(which="minor", bottom=False, left=False)

    return im, cbar

def annotate_heatmap(im, data=None, valfmt="{x:.1f}",
                     textcolors=["black", "white"],
                     threshold=None, **textkw):
    """
    A function to annotate a heatmap.

    Arguments:
        im         : The AxesImage to be labeled.
    Optional arguments:
        data       : Data used to annotate. If None, the image's data is used.
        valfmt     : The format of the annotations inside the heatmap.
                     This should either use the string format method, e.g.
                     "$ {x:.2f}", or be a :class:`matplotlib.ticker.Formatter`.
        textcolors : A list or array of two color specifications. The first is
                     used for values below a threshold, the second for those
                     above.
        threshold  : Value in data units according to which the colors from
                     textcolors are applied. If None (the default) uses the
                     middle of the colormap as separation.

    Further arguments are passed on to the created text labels.
    """

    if not isinstance(data, (list, np.ndarray)):
        data = im.get_array()

    # Normalize the threshold to the images color range.
    if threshold is not None:
        threshold = im.norm(threshold)
    else:
        threshold = im.norm(data.max())/2.

    # Set default alignment to center, but allow it to be
    # overwritten by textkw.
    kw = dict(horizontalalignment="center",
              verticalalignment="center")
    kw.update(textkw)

    # Get the formatter in case a string is supplied
    if isinstance(valfmt, str):
        valfmt = matplotlib.ticker.StrMethodFormatter(valfmt)

    # Loop over the data and create a `Text` for each "pixel".
    # Change the text's color depending on the data.
    texts = []
    for i in range(data.shape[0]):
        for j in range(data.shape[1]):
            kw.update(color=textcolors[im.norm(data[i, j]) > threshold])
            text = im.axes.text(j, i, "{0:.1f}".format(data[i, j]), **kw)
            texts.append(text)

    return texts


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

import numpy as np
import matplotlib.pyplot as plt

import pycity_base.classes.supply.BES as BES
import pycity_base.classes.supply.Boiler as Boiler

import pycity_base.classes.supply.CHP as CHP
import pycity_calc.economic.energy_sys_cost.chp_cost as chp_cost
import pycity_calc.economic.energy_sys_cost.tes_cost as tes_cost
import pycity_calc.economic.energy_sys_cost.boiler_cost as boiler_cost

import pycity_calc.energysystems.heatPumpSimple as HP

def dim_decentral_hp(environment, sh_curve, t_biv=-2, hp_type='aw', lowerActivationLimit=0.5, tSink=35, temp_norm=-12, t_blocked=4):
    """
    Dimensioning of air/water and brine/water heat pump.

    Parameters
    ----------
    environment : Environment class of city-object
    sh_curve : Space heating curve of building
    t_biv : bivalence temperature
    hp_type : type of hp, 'aw'/'ww'
    lowerActivationLimit : lower activation limit of heat pump
    tSink : minimal flow temperature
    temp_norm : -12°C for Aachen

    Returns
    -------
    q_hp_biv : float
        Power of heat pump in bivalence point (-> Q_nom)
    tMax : int
        max. flow temperature
    lowerActivationLimit : float
        lower activation limit of heat pump
    tSink : int
        return flow temperature
    t_dem_ldc : list
        space heating demand sorted by outside temperature
    biv_ind : int
        index of bivalence point in t_dem_ldc
    """

    if hp_type == 'aw':
        # monoenergetic operation mode

        # get ldc with sh-demand and tAmbient
        t_dem_ldc = get_t_demand_list(temp_curve=environment.weather.tAmbient,
                                      th_curve=sh_curve)

        # plot temp-sorted demand curve
        # plt.plot(sorted(environment.weather.tAmbient), t_dem_ldc)

        # Find bivalence point and necessary thermal power of heat pump
        for i in range(len(t_dem_ldc)):

            # iterate through temp from min to max until temp is bigger than t_biv
            if sorted(environment.weather.tAmbient)[i] > t_biv:

                # get index of bivalence point in sorted temp curve
                biv_ind = i - 1

                # th. power in bivalence point (straight line between max(sh-demand) at temp_norm and 0 kW at 20°C)
                q_hp_biv = max(sh_curve) / (temp_norm - 20) * (t_biv - 20)

                # calculate produced heat by peak supply
                q2 = 0
                for q in sh_curve:
                    if q > q_hp_biv:
                        q2 += q - q_hp_biv

                # If curves are generated with method 1 (SLP by Hellwig), ratio is always 0.91% - maybe bc curves are only scaled up?
                print('peak supply produces ' + str(np.round(q2*100/np.sum(sh_curve),2)) + '% of annual heat demand.')

                ee_ratio = 1 - q2 / np.sum(sh_curve)

                break
        else:
            raise Exception('Error in demand calculation at bivalence point')

        # Choose heat pump from catalogue
        q_nom, cop_list, tMax, tSink = choose_hp(q_ideal=q_hp_biv,
                                                 t_biv=t_biv)

        return q_nom, cop_list, tMax, lowerActivationLimit, tSink, t_dem_ldc, biv_ind, ee_ratio

    elif hp_type == 'ww':
        # monovalent operation mode

        # two compressors
        lowerActivationLimit = 0.5

        # factor for blocked time (energy supply companies)
        f = round(24/(24-t_blocked), 1)

        # calculate necessary th. power
        q_hp = max(sh_curve) * f

        # choose heat pump from catalogue
        q_nom, cop_list, tMax, tSink = choose_hp(q_ideal=q_hp,
                                                 method=0,
                                                 hp_type='ww',
                                                 t_ground=environment.temp_ground)

        return q_nom, cop_list, tMax, lowerActivationLimit, tSink

    else:
        raise Warning('wrong heatpump type defined.')

def choose_hp(q_ideal, t_biv=-2, method=0, hp_type='aw', t_ground=10):
    """
    Choose heat pump depending on desired thermal power and bivalence point
    Method 0: Return best possible Heat Pump with average COPs
    Method 1: Choose heat pump from catalogue of existing devices (source: Dimplex)

    Parameters
    ----------
    q_ideal
    t_biv
    method
    hp_type
    t_ground

    Returns
    -------

    """

    # choose heat pump with ideal th.power and average cop
    if method == 0:

        # Air/Water heatpump
        if hp_type == 'aw':
            best_q = q_ideal

            # COPs for [A-7/W35, A2/W35, A7/W35]
            best_cop_list = [2.9, 3.7, 4.4]

            # 55/35 normal temp spread (Vor-/Ruecklauf)
            tMax = 43   #55
            tSink = 43  #35

        # Brine/Water heatpump
        elif hp_type == 'ww':
            best_q = q_ideal

            # COPs depending on constant t_ground (depending on depth of collector)
            if t_ground == 0:
                best_cop_list = [4.8]
            elif t_ground == 10:
                best_cop_list = [6]
            else:
                raise Warning('ground temperature not supported: ' + str(t_ground))

            tMax = 62
            tSink = 35

        else:
            raise Warning('heat pump type unknown: ' + hp_type)

        return best_q, best_cop_list, tMax, tSink

    # choose heatpump from catalogue
    elif method == 1:

        # Air/Water heatpump
        if hp_type == 'aw':

            # [heating power],[COP] (at A-7/W35, A2/W35, A7/W35), [tMax, tMin] (source: Dimplex)
            hp_list = {'LA6TU':([4000,5100,6400],[2.9,3.8,4.6],[60,18]),'LA9TU':([5200,7500,9200],[2.8,3.6,4.2],[58,18]),
                       'LA12TU':([7600,9400,11600],[2.9,3.7,4.3],[58,18]),'LA17TU':([10300,14600,19600],[2.9,3.7,4.4],[58,18]),
                       'LA25TU':([16700,19600,26100],[3.0,3.7,4.4],[58,18]),'LA40TU':([23800,30000,35700], [3.0,3.8,4.4],[58,18])}

            # get data for smallest heatpump in catalogue
            best_hp = 'LA6TU'
            best_q_biv = (hp_list['LA6TU'][0][2] - hp_list['LA6TU'][0][0]) / (7 + 7) * (t_biv + 7) + hp_list['LA6TU'][0][0]
            best_cop_list = hp_list['LA6TU'][1]

            # iterate through catalogue and choose most suitable heatpump
            for hp in hp_list:

                # get th.power list of current hp
                q_nom = hp_list[hp][0]

                # get th.power at t_biv for current hp (interpolation)
                q_biv = (q_nom[2] - q_nom[0])/(7 + 7) * (t_biv + 7) + q_nom[0]

                # update best_hp if current hp is more suitable
                if abs(q_biv - q_ideal) <= abs(best_q_biv - q_ideal):
                    best_hp = hp
                    best_q_biv = q_biv
                    best_cop_list = hp_list[hp][1]

            print('Best A/W-HP: ' + best_hp)

            # temperature constraints
            tMax = hp_list[best_hp][2][0]
            tSink = 35

            return best_q_biv, best_cop_list, tMax, tSink

        elif hp_type == 'ww':
            # [heating power],[COP] (at B0/W35, B10/W35),[tMax] (source: Dimplex)
            hp_list = {'SI8TU':([8100,11000],[4.8,6.2],[62]),'SI11TU':([10900,14100],[4.9,6.2],[62]),'SI14TU':([13900,18000],[5,6.2],[62]),
                       'SI18TU':([17500,22000],[4.7,5.6],[62]),'SI22TU':([22900,28000],[4.4,5.2],[58])}

            if t_ground == 10:
                t_ind = 0
            elif t_ground == 10:
                t_ind = 1
            else:
                raise Warning('hp for specific ground temperature not available: ' + str(t_ground) + '°C')

            # iterate through catalogue and choose most suitable heatpump
            for hp in hp_list:
                q_nom = hp_list[hp][0][t_ind]

                # choose first hp with th.power higher than needed
                if q_nom >= q_ideal:
                    best_hp = hp
                    best_q = q_nom
                    best_cop_list = hp_list[hp][1]
                    break
            else:
                raise Warning('Heat pumps in catalogue too small. Needed ' + str(round(q_ideal/1000,2))+'kW'+'; Found: '
                              + str(round(hp_list['SI22TU'][0][t_ind]/1000,2)) + ' kW.')

            print('Best S/W-HP: ' + best_hp)

            # temperature constraints
            tMax = hp_list[best_hp][2][0]
            tSink = 35

            return best_q, best_cop_list, tMax, tSink

        else:
            raise Warning('heat pump type unknown: ' + hp_type)

    else:
        raise Warning('method unknown: ' + str(method))

def dim_decentral_chp(th_LDC, q_total, method=0):
    """
    Dimensioning of CHP in decentralized supply system.

    Parameters
    ----------
    th_LDC : list : thermal load demand curve
    q_total : float : total heat demand per year in Wh
    method : int
        0 : method according to krimmling 2011
        1 : method to comply with EEWaermeG

    Returns
    -------
    eta_el : el. efficiency
    eta_th : th. efficiency
    p_nom : el. nominal power
    q_nom : th. nominal power
    t_x : full load hours
    t_ann_op : annual operation hours
    """

    # ---------------- Method 0: Krimmling and interviews ----------------
    if method == 0:

        # Set minimum th. power for chp
        q_chp = 0.1 * max(th_LDC)

        # Choose CHP from catalogue with minimum th. power and get estimated op.time (annual op. hours, full load hours)
        [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)
        (t_ann_op_max, t_x_max) = get_chp_ann_op_time(q_nom, th_LDC)

        # Set minimum full load hours
        chp_flh = 5000

        # get th. power with minimal full load hours (max th. power for chp)
        q_chp = th_LDC[chp_flh]

        # Choose CHP from catalogue with minimum th. power and get estimated op.time
        [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)
        (t_ann_op_min, t_x_min) = get_chp_ann_op_time(q_nom, th_LDC)

        # search best solution between chp with min and max full load hours
        if t_x_min <= t_x_max:

            # dummy values
            best_costs = 1000000000
            best_sol_costs = 0, 0, 0, 0, 0, 0

            array_q_nom = []
            array_t_x = []
            array_t_ann_op = []
            array_deckung = []
            array_costs = []

            # search for most economical chp - each iteration increases flh by 20h
            for t in range(t_x_min, t_x_max+1, 20):  # t_x_max+1 to have at least 1 value in range if t_max = t_min
                cost_cap = []
                cost_dem = []
                rev = []

                # choose best chp at timestep t
                q_chp = th_LDC[t]
                [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)

                # get estimated operational hours
                (t_ann_op, t_x) = get_chp_ann_op_time(q_nom, th_LDC)
                chp_ratio = q_nom * t_ann_op / q_total

                # Calculate estimated costs for chp
                q_chp_nom = q_nom / 1000  # in kW
                p_chp_nom = p_nom / 1000  # in kW

                a_chp = 0.08 * (1 + 0.08) ** 15 / ((1 + 0.08) ** 15 - 1)  # annuity factor
                chp_invest = chp_cost.calc_invest_cost_chp(p_el_nom=p_chp_nom,
                                                           method='asue2015',
                                                           with_inst=True,
                                                           use_el_input=True,
                                                           q_th_nom=None)
                cost_cap.append(chp_invest * a_chp)
                cost_dem.append((q_chp_nom) / eta_th * t_ann_op * 0.0661)
                rev.append(p_chp_nom * t_ann_op * get_el_feedin_tariff_chp(q_chp_nom))

                # Calculate estimated costs for thermal energy storage
                v_tes = q_chp_nom * 60  # in liter

                a_tes = 0.08 * (1 + 0.08) ** 20 / ((1 + 0.08) ** 20 - 1)
                tes_invest = tes_cost.calc_invest_cost_tes(v_tes / 1000, method='spieker')

                cost_cap.append(tes_invest * a_tes)
                rev.append(get_subs_tes_chp(chp_ratio, v_tes, tes_invest, p_chp_nom))
                rev.append(get_subs_minichp(p_chp_nom, q_chp_nom, v_tes))

                # Calculate costs for Boiler
                q_boiler_nom = max(max(th_LDC) - q_nom, 0) / 1000

                a_boiler = 0.08 * (1 + 0.08) ** 18 / ((1 + 0.08) ** 18 - 1)
                boiler_invest = boiler_cost.calc_abs_boiler_cost(q_boiler_nom, method='viess2013')

                cost_cap.append(boiler_invest * a_boiler)
                cost_dem.append((q_total / 1000 - q_chp_nom * t_ann_op) / 0.8 * 0.0661)

                costs_total = sum(cost_cap) + sum(cost_dem) - sum(rev)

                array_q_nom.append(q_nom)
                array_t_x.append(t_x)
                array_t_ann_op.append(t_ann_op)
                array_deckung.append((t_ann_op * q_nom) / q_total)
                array_costs.append(costs_total)

                # Choose current option if best option so far
                if costs_total < best_costs:
                    best_sol_costs = (eta_el, eta_th, p_nom, q_nom, t_x, t_ann_op)
                    best_costs = costs_total

            (eta_el, eta_th, p_nom, q_nom, t_x, t_ann_op) = best_sol_costs

        else:
            return None

    # ---------------- Method 1: EEWaermeG ----------------

    elif method == 1:

        eewg = False
        count = 0

        # Requirement are at least 50% heat per year from chp
        # Calculate sufficient th.power for CHP for average operation
        q_chp = 0.5 * q_total / 8760

        # Choose appropriate device and calculate estimated operation time
        [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)
        (t_ann_op, t_x) = get_chp_ann_op_time(q_nom, th_LDC)

        while not eewg:
            # Calculate ratio of chp heat usage (min. 50% necessary)
            ee_ratio = q_nom * t_ann_op / q_total

            if ee_ratio > 0.5:

                # calculate primary energy for separate production
                refeta_th = 0.85  # th. reference efficiency rate for devices older than 2016 (natural gas)
                refeta_el = 0.525  # el. reference efficiency rate for devices built between 2012 and 2015 (natural gas)

                # Primary energy savings in % calculated according to EU directive 2012/27/EU
                pee = (1 - 1 / ((eta_th / refeta_th) + (eta_el / refeta_el))) * 100

                # PE savings for CHP with el.power >= 1 MW have to be at least 10%
                if p_nom >= 1000000:
                    if pee >= 10:
                        eewg = True
                        break

                # PE savings for CHP with el. power < 1 MW have to be existent
                else:
                    if pee > 0:
                        eewg = True
                        break

                if ee_ratio >= 1:
                    raise Warning('EEWaermeG not satisfied: unrealistic values! (Q_chp >= Q_total)')

            # Raise th.power with 5% per iteration if conditions for EEWaermeG were not met
            q_chp = 0.5 * q_total / 8760 + count * q_total / (8760 * 100)
            [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)
            (t_ann_op, t_x) = get_chp_ann_op_time(q_nom, th_LDC)
            count += 5
    else:
        raise Exception('Invalid input for method!')

    return eta_el, eta_th, p_nom, q_nom, t_x, t_ann_op

def dim_central_chp(th_LDC, q_total, method=0):
    """
    Dimensioning of CHP in centralized supply system.

    Parameters
    ----------
    method : int
        0 : Standard approach (Krimmling 2011):
                - Q_th > 10% (< 20%) of max. thermal demand
                - Full-Load-hours preferably > 6000 h/a (minimum 5000 h/a)
        1 : Approach to satisfy EEWaermeG without any other constraints
        2 : Choose device where BAFA-Support for heating networks is available (Q_chp >= 75% of total annual heat
            demand) as described in KWKG without any other constraints

    Returns
    -------
    chp : ChpExtended
    """

    # ---------------- Method 0: Krimmling and interviews ----------------
    if method == 0:

        # Set minimum th. power for chp
        q_chp = 0.1*max(th_LDC)

        # Choose CHP from catalogue with minimum th. power and get estimated op.time (annual op. hours, full load hours)
        [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)
        (t_ann_op_max, t_x_max) = get_chp_ann_op_time(q_nom, th_LDC)

        # Set minimum full load hours
        chp_flh = 5000

        # get th. power with minimal full load hours (max th. power for chp)
        q_chp = th_LDC[chp_flh]

        # Choose CHP from catalogue with minimum th. power and get estimated op.time
        [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)
        (t_ann_op_min, t_x_min) = get_chp_ann_op_time(q_nom, th_LDC)

        # search best solution between chp with min and max full load hours
        if t_x_min < t_x_max:

            # dummy values
            best_costs = 1000000000
            best_sol_costs = 0,0,0,0,0,0

            array_q_nom = []
            array_t_x = []
            array_t_ann_op = []
            array_deckung = []
            array_costs = []

            # search for most economical chp - each iteration increases flh by 20h
            for t in range(t_x_min, t_x_max+1, 20):  # t_x_max+1 to have at least 1 value in range if t_max = t_min
                cost_cap = []
                cost_dem = []
                rev = []

                # choose best chp at timestep t
                q_chp = th_LDC[t]
                [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)

                # get estimated operational hours
                (t_ann_op, t_x) = get_chp_ann_op_time(q_nom, th_LDC)
                chp_ratio = q_nom * t_ann_op / q_total

                # Calculate estimated costs for chp
                q_chp_nom = q_nom/1000  # in kW
                p_chp_nom = p_nom/1000  # in kW

                a_chp = 0.08 * (1 + 0.08) ** 15 / ((1 + 0.08) ** 15 - 1)    # annuity factor
                chp_invest = chp_cost.calc_invest_cost_chp(p_el_nom=p_chp_nom,
                                                           method='asue2015',
                                                           with_inst=True,
                                                           use_el_input=True,
                                                           q_th_nom=None)
                cost_cap.append(chp_invest * a_chp)
                cost_dem.append((q_chp_nom)/eta_th * t_ann_op * 0.0661)
                rev.append(p_chp_nom * t_ann_op * get_el_feedin_tariff_chp(q_chp_nom))

                # Calculate estimated costs for thermal energy storage
                v_tes = q_chp_nom * 60  # in liter

                a_tes = 0.08 * (1 + 0.08) ** 20 / ((1 + 0.08) ** 20 - 1)    # annuity factor
                tes_invest = tes_cost.calc_invest_cost_tes(v_tes/1000, method='spieker')

                cost_cap.append(tes_invest * a_tes)
                rev.append(get_subs_tes_chp(chp_ratio, v_tes, tes_invest, p_chp_nom))
                rev.append(get_subs_minichp(p_chp_nom, q_chp_nom, v_tes))

                # Calculate costs for Boiler
                q_boiler_nom = max(max(th_LDC) - q_nom, 0)/1000

                a_boiler = 0.08 * (1 + 0.08) ** 18 / ((1 + 0.08) ** 18 - 1)    # annuity factor
                boiler_invest = boiler_cost.calc_abs_boiler_cost(q_boiler_nom, method='viess2013')

                cost_cap.append(boiler_invest * a_boiler)
                cost_dem.append((q_total/1000 - q_chp_nom * t_ann_op)/0.8*0.0661)  # operational costs for boiler

                costs_total = sum(cost_cap) + sum(cost_dem) - sum(rev)

                array_q_nom.append(q_nom)
                array_t_x.append(t_x)
                array_t_ann_op.append(t_ann_op)
                array_deckung.append((t_ann_op * q_nom) / q_total)
                array_costs.append(costs_total)

                # Choose current option if best option so far
                if costs_total < best_costs:
                    best_sol_costs = (eta_el, eta_th, p_nom, q_nom, t_x, t_ann_op)
                    best_costs = costs_total

            (eta_el, eta_th, p_nom, q_nom, t_x, t_ann_op) = best_sol_costs

        else:
            return None

    # ---------------- Method 1: EEWaermeG ----------------

    elif method == 1:

        eewg = False
        count = 0

        # Requirement are at least 50% heat per year from chp
        # Calculate sufficient th.power for CHP for average operation
        q_chp = 0.5 * q_total / 8760

        # Choose appropriate device and calculate estimated operation time
        [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)
        (t_ann_op, t_x) = get_chp_ann_op_time(q_nom, th_LDC)

        while not eewg:
            # Calculate ratio of chp heat usage (min. 50% necessary)
            ee_ratio = q_nom * t_ann_op / q_total

            if ee_ratio > 0.5:

                # calculate primary energy for separate production
                refeta_th = 0.85  # th. reference efficiency rate for devices older than 2016 (natural gas)
                refeta_el = 0.525  # el. reference efficiency rate for devices built between 2012 and 2015 (natural gas)

                # Primary energy savings in % calculated according to EU directive 2012/27/EU
                pee = (1 - 1 / ((eta_th / refeta_th) + (eta_el / refeta_el))) * 100

                # PE savings for CHP with el.power >= 1 MW have to be at least 10%
                if p_nom >= 1000000:
                    if pee >= 10:
                        eewg = True
                        break

                # PE savings for CHP with el. power < 1 MW have to be existent
                else:
                    if pee > 0:
                        eewg = True
                        break

                if ee_ratio >= 1:
                    print('EEWaermeG not satisfied: unrealistic values! (Q_chp >= Q_total)')
                    break

            # Raise th.power with 5% per iteration if conditions for EEWaermeG were not met
            q_chp = 0.5 * q_total / 8760 + count * q_total / (8760 * 100)
            [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)
            (t_ann_op, t_x) = get_chp_ann_op_time(q_nom, th_LDC)
            count += 5

    # ---------------- Method 2: BAFA support ----------------

    elif method == 2:

        # best possible full-load-hours (top border)
        chp_flh = 7000

        # Calculate CHP th.power for max full load hours
        q_chp = th_LDC[chp_flh]

        # Choose appropriate device and calculate estimated operation time
        [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)
        (t_ann_op, t_x) = get_chp_ann_op_time(q_nom, th_LDC)

        while True:
            # BAFA subsidy applicable if 75% of annual heat demand are supplied by CHP
            if (q_nom * t_ann_op / q_total) >= 0.75:

                print('CHP: BAFA funding for LHN possible! Coverage: ' + str(
                    round(q_nom * t_ann_op * 100 / q_total, 2)) + '%')
                print('Full load hours:', t_x, 'h/a; Annual operation time:', t_ann_op, 'h/a')
                break

            # if annual CHP heat supply ratio is smaller than 75%, lower flh (raises CHP th.power)
            else:
                chp_flh -= 75

                if chp_flh > 1000:
                    q_chp = th_LDC[chp_flh]
                    [eta_el, eta_th, p_nom, q_nom] = choose_chp(q_chp)
                    (t_ann_op, t_x) = get_chp_ann_op_time(q_nom, th_LDC)
                    print(round(q_nom * t_ann_op * 100 / q_total,2),'%')

                # break if full load hours are below 1000 h/yr
                else:
                    print('CHP: BAFA funding for LHN not possible!')
                    return None

    # unknown method
    else:
        raise Exception('Invalid input for method!')

    return eta_el, eta_th, p_nom, q_nom, t_x, t_ann_op

def get_LDC(curve):
    """
    Calculates load duration curve (Jahresdauerlinie)

    Parameters
    ----------
    curve : thermal or electrical curve

    Returns
    -------
    sorted curve : sorted load duration curve
    """
    return sorted(curve, reverse=True)

def get_t_demand_list(temp_curve, th_curve):
    """
    Sorts thermal energy demand based on values of ambient temperatures.

    Parameters
    ----------
    temp_curve : list of ambient temperatures for one year
    th_curve : list of thermal energy demand for one year

    Returns
    -------
    t_demand_curve : thermal energy curve sorted based on ambient temperature values
    """
    return [th_demand for _, th_demand in sorted(zip(temp_curve, th_curve))]

def get_chp_ann_op_time(q_nom, th_LDC):
    """
    Calculates annual operation time and full load hours of chp

    Parameters
    ----------
    q_nom   : Nominal thermal power of chp
    th_LDC  : thermal Load duration curve

    Returns
    -------
    t_ann_operation : total annual operation time of chp
    t_x             : Full load hours of chp
    """

    # calculate annual operation time according to Krimmling 2011
    for q_m in th_LDC:

        # find crossing point on LDC (indicates full load hours)
        if q_m <= q_nom:
            t_x = th_LDC.index(q_m)
            break
    else:
        print('Warning! No crossing point between CHP th.power and LDC found. Check parameters or add devices!')
        return 0, 0

    # Find point in time, when area between t:t_x and t_x:8760 are the same size (see Krimmling p.131, Abb.4-15)

    # difference between areas
    delta_a = 8760 * th_LDC[0]

    # make ldc numpy array
    ldc = np.array(th_LDC)

    # iterate through timesteps between full load hours of chp and end of year
    for t in range(t_x, 8760):

        # area 1
        a1 = q_m * (t - t_x) - np.sum(ldc[t_x:t])
        # area 2
        a2 = np.sum(ldc[t:8760])

        # if difference between areas is bigger than last iteration
        if delta_a <= abs(a2 - a1):

            # t-1 indicates index where a1 is (almost) a2 -> estimated annual hours of operation
            t_ann_operation = t - 1
            return t_ann_operation, t_x

        else:
            delta_a = a2 - a1

def choose_chp(q_ideal, method=0):
    """
    Choose CHP device
    Method 0: Return device with desired thermal power and average effiency
    Method 1: Out of catalogue of available devices

    Parameters
    ----------
    q_ideal : exact thermal power needed of chp

    Returns
    -------
    specs : list with specifications of chp device [eta_el, eta_th, p_nom, q_nom]
    """

    if method == 0:
        q_nom = q_ideal
        #eta_el = 0.27  # average for real devices
        #eta_th = 0.65  # average for real devices
        eta_el = 0.45   # default in pyCity_calc
        eta_th = 0.45   # default in pyCity_calc

        # Calculate el.power based on efficiency and ideal th.power
        p_nom = round(q_nom/eta_th * eta_el/10)*10

        return eta_el, eta_th, p_nom, q_nom

    else:
        # TODO: Import data from excel
        # [eta_el, eta_th, p_nom, q_nom] - (Source: BHKW-Kenndaten 2014, S.26)
        chp_list = {'vai1': [0.263, 0.658, 1000, 2500], 'intelli': [0.245, 0.613, 2600, 6500],
                    'vai2': [0.25, 0.667, 3000, 8000], 'kirsch': [0.19, 0.76, 1900, 9000],
                    'vai3': [0.247, 0.658, 4700, 12500], 'ecp': [0.286, 0.644, 6000, 13500],
                    'vie': [0.27, 0.671, 6000, 14900], 'hoef': [0.27, 0.664, 7000, 17200],
                    'rmb7.2': [0.263, 0.657, 7200, 18000],'oet8': [0.268, 0.633, 8000, 19000],
                    'xrgi9': [0.289, 0.641, 9000, 20000], 'kwe': [0.268, 0.818, 7500, 22900],
                    'rmb11.0': [0.289, 0.632, 11000, 24000], 'oet12': [0.279, 0.605, 12000, 26000],
                    'motat': [0.279, 0.651, 12000, 28000],'xrgi15': [0.307, 0.613, 15000, 30000],
                    'ews': [0.711, 0.311, 14000, 32000], 'asv1534': [0.306, 0.694, 15000, 34000],
                    'sb16': [0.314, 0.72, 16000, 36700], 'enrtc': [0.323, 0.629, 20000, 39000],
                    'xrgi20': [0.32, 0.64, 20000, 40000]}

        # define starting specs
        if q_ideal < 2500:
            specs = chp_list['vai1']
        else:
            specs = [0, 0, 0, 0]

        # choose most suitable chp regarding th.power
        for dev in chp_list.values():
            q_nom = dev[3]
            if abs(q_nom-q_ideal) < abs(specs[3]-q_ideal):
                specs = dev[:]

        return specs

def get_el_feedin_tariff_chp(q_nom, el_feedin_epex=0.02978, vnn=0.01):
    """
    Calculates feed-in tariff for CHP-produced electricity.

    Parameters
    ----------
    q_nom : nominal thermal power of chp in kW
    el_feedin_epex : epex price for electricity in Euro/kWh
    vnn : avoided grid charges ("vermiedenes Netznutzungsentgelt") in Euro/kWh

    Returns
    -------
    feed-in tariff in EUR/kWh
    """

    # KWKG 2016 revenues for el. feed-in + feedin price from epex + avoided grid charges
    if q_nom < 50:
        return 0.08+el_feedin_epex+vnn  # Euro/kWh, only paid for 60.000 flh
    elif 50 <= q_nom < 100:
        return 0.06+el_feedin_epex+vnn  # Euro/kWh, only paid for 30.000 flh
    elif 100 <= q_nom < 250:
        return 0.05+el_feedin_epex+vnn  # Euro/kWh, only paid for 30.000 flh
    elif 250 <= q_nom < 2000:
        return 0.044+el_feedin_epex+vnn  # Euro/kWh, only paid for 30.000 flh
    else:  # q_nom > 2000:
        return 0.031+el_feedin_epex+vnn  # Euro/kWh, only paid for 30.000 flh

def get_subs_minichp(p_nom, q_nom, v_tes):
    """
    Calculate BAFA subsidy for mini chp devices.

    Parameters
    ----------
    p_nom : chp el. power in kW
    q_nom : chp th. power in kW
    v_tes : volume of thermal energy storage in liter

    Returns
    -------
    bafa_subs_chp : subsidy for mini-chp
    """

    # BAFA subsidy for Mini-CHP (CHP device must be listed on BAFA-list)
    bafa_subs_chp = 0

    if v_tes > 0:

        # el.power smaller than 20 kW and TES-volume at least 60 l/kW_th
        if p_nom < 20 and v_tes/q_nom >= 60:

            if p_nom < 1:
                bafa_subs_chp = 1900

            elif p_nom < 4:
                bafa_subs_chp = 1900 + (p_nom / 1000 - 1) * 300

            elif p_nom < 10:
                bafa_subs_chp = 1900 + 3 * 300 + (p_nom / 1000 - 4) * 100

            else:  # bes.chp.pNominal < 20:
                bafa_subs_chp = 1900 + 3 * 300 + 6 * 100 + (p_nom / 1000 - 10) * 10

    return bafa_subs_chp

def get_subs_tes_chp(chp_ratio, v_tes, tes_invest, p_nom):
    """
    Calculate KWKG subsidy for TES in combination with CHP

    Parameters
    ----------
    chp_ratio : chp_heat/total_heat per year
    v_tes : tes volume in liter
    tes_invest
    p_nom : el. power in kW

    Returns
    -------
    kwkg_subs_tes : amount of subsidy
    """

    kwkg_subs_tes = 0

    # TES volume in m3
    v_m3 = v_tes/1000

    # subsidy applies only for systems with ratio of chp heat per year higher than 50%
    if chp_ratio >= 0.5:

        if 1.0 <= v_m3 <= 50.0:
            kwkg_subs_tes = 250 * v_m3

        elif v_m3 > 50:
            kwkg_subs_tes = 0.3 * tes_invest

        else:
            if v_m3 >= 0.3 * p_nom:
                kwkg_subs_tes = 250 * v_m3

    return kwkg_subs_tes

def get_bafa_subs_hp(q_nom, spf):
    """
    BAFA subsidy for a/w heat pumps with nominal power <= 100 kW and spf >= 3.5

    heat pump must be used for:
        - combined space heating and dhw
        - only space heating if dhw is generated with renewable energy
        - generating heat for lhn

    Values only for retrofit; subsidy for new buildings is lower (not implemented)

    Parameters
    ----------
    q_nom : thermal nominal power of heat pump in kW
    spf : seasonal performance factor ("Jahresarbeitszahl")

    Returns
    -------
    subs : amount of subsidy
    """
    subs = 0

    # innovation subsidy
    if q_nom <= 100 and spf >= 4.5:
        if q_nom < 32.5:
            subs = 1950
        else:
            subs = 60*q_nom

    # base subsidy
    if q_nom <= 100 and spf >= 3.5:
        if q_nom < 32.5:
            subs = 1300
        else:
            subs = 40*q_nom

    return subs

def get_kfw_subs_430(invest_cost):
    '''
    KfW subsidy for single measures regarding program 430. Currently only applicable for first connection to lhn.
    For further constraints see KfW website.

    Parameters
    ----------
    invest_cost : invest cost of single measure

    Returns
    -------
    kfw_subs_430 : possible subsidy from KfW
    '''

    kfw_subs_430 = invest_cost * 0.1

    return kfw_subs_430


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

Drive the landlab potentiality flow routing component.

Created on Wed Mar 4 2015

@author: danhobley
"""
from __future__ import print_function

import numpy as np
from pylab import close, figure, show

from landlab import ModelParameterDictionary, RasterModelGrid
from landlab.components import FastscapeEroder, FlowAccumulator, PotentialityFlowRouter
from landlab.plot.imshow import imshow_grid

nrows = 100
ncols = 100
dx = 1000.

close("all")
mg = RasterModelGrid((nrows, ncols), dx)
z = (3000. - mg.node_x) * 0.5
# z = -mg.node_x-mg.node_y
# z = np.sqrt(mg.node_x**2 + mg.node_y**2)
# z = mg.node_x + mg.node_y
# val_to_replace_with = z.reshape((nrows,ncols))[3,3]
# z.reshape((nrows,ncols))[2:5,:] = val_to_replace_with
mg.at_node["topographic__elevation"] = z

# mg.set_fixed_value_boundaries_at_grid_edges(True, True, True, True)
mg.set_closed_boundaries_at_grid_edges(False, True, True, True)
figure(3)
imshow_grid(mg, mg.status_at_node)

mg.at_node["water__unit_flux_in"] = np.ones_like(z)

pfr = PotentialityFlowRouter(mg, "pot_fr_params.txt")

pfr.route_flow(route_on_diagonals=True)

figure(1)
imshow_grid(mg, "surface_water__discharge")
figure(2)
imshow_grid(mg, "topographic__elevation")

out_sum = np.sum(mg.at_node["surface_water__discharge"].reshape((nrows, ncols))[-3, :])
print(out_sum)
print(
    np.sum(mg.at_node["water__unit_flux_in"]),
    np.sum(mg.at_node["surface_water__discharge"][mg.boundary_nodes]),
)
print(out_sum - np.sum(mg.at_node["water__unit_flux_in"]))

show()

print(mg.at_node["surface_water__discharge"].reshape((nrows, ncols)))
print(
    np.sum(mg.at_node["surface_water__discharge"].reshape((nrows, ncols)), axis=0) / 3.
)

# now a run with a grid...

# mg = RasterModelGrid((nrows, ncols), dx)
# # make a topo to test on:
# mg.at_node['topographic__elevation'] = np.zeros(mg.number_of_nodes)
# fr = FlowRouter(mg)
# fsp = FastscapeEroder(mg, './pot_fr_params.txt')
# inputs = ModelParameterDictionary('./pot_fr_params.txt')
# dt = inputs.read_float('dt')
# time_to_run = inputs.read_float('run_time')
# uplift = inputs.read_float('uplift_rate')
# mg.set_fixed_value_boundaries_at_grid_edges(True, True, True, True)
# # make some K values in a field to test
# mg.at_node['K_values'] = 0.1+np.random.rand(nrows*ncols)/10.
# # elapsed_time = 0. #total time in simulation
# # while elapsed_time < time_to_run:
# #     #add uplift
# #     mg.at_node['topographic__elevation'][mg.core_nodes] += uplift*dt
# #     print elapsed_time
# #     if elapsed_time+dt>time_to_run:
# #         print "Short step!"
# #         dt = time_to_run - elapsed_time
# #     mg = fr.route_flow(grid=mg)
# #     mg,_,_ = sp.erode(mg)
# #     elapsed_time += dt

inputs = ModelParameterDictionary("./pot_fr_params.txt")
nrows = 200  # inputs.read_int('nrows')
ncols = 200  # inputs.read_int('ncols')
dx = inputs.read_float("dx")
dt = inputs.read_float("dt")
time_to_run = inputs.read_float("run_time")
# nt needs defining
uplift = inputs.read_float("uplift_rate")
init_elev = inputs.read_float("init_elev")

mg = RasterModelGrid((nrows, ncols), dx)

# create the fields in the grid
mg.add_zeros("topographic__elevation", at="node")
z = mg.zeros(at="node") + init_elev
mg["node"]["topographic__elevation"] = z + np.random.rand(len(z)) / 1000.

# make some K values in a field to test
mg.at_node["K_values"] = 0.00001 + np.random.rand(nrows * ncols) / 100000.

# mg.at_node['water__unit_flux_in'] = dx*dx*np.ones_like(z)
mg.at_node["water__unit_flux_in"] = (
    dx * dx * np.ones_like(z) * 100. / (60. * 60. * 24. * 365.25)
)  # remember, flux is /sec, so this is a small number!
# mg.set_closed_boundaries_at_grid_edges(False, False, True, True)
# mg.set_closed_boundaries_at_grid_edges(True, False, True, True)

print("Running ...")

# instantiate the components:
fr = FlowAccumulator(mg, flow_director="D8")
# load the Fastscape module too, to allow direct comparison
fsp = FastscapeEroder(mg, "./pot_fr_params.txt")

# perform the loop:
elapsed_time = 0.  # total time in simulation
while elapsed_time < time_to_run:
    print(elapsed_time)
    if elapsed_time + dt > time_to_run:
        print("Short step!")
        dt = time_to_run - elapsed_time
    mg = fr.run_one_step()
    # print 'Area: ', numpy.max(mg.at_node['drainage_area'])
    # mg = fsp.erode(mg)
    mg = fsp.erode(mg, K_if_used="K_values")
    # mg,_,_ = sp.erode(mg, dt, node_drainage_areas='drainage_area', slopes_at_nodes='topographic__steepest_slope')
    # add uplift
    mg.at_node["topographic__elevation"][mg.core_nodes] += uplift * dt
    elapsed_time += dt

pfr = PotentialityFlowRouter(mg, "pot_fr_params.txt")
pfr.route_flow(return_components=True)  # route_on_diagonals=False)

figure("Topo")
imshow_grid(mg, "topographic__elevation")
figure("Potentiality flow fluxes")
imshow_grid(mg, "surface_water__discharge")
figure("D8 drainage areas")
imshow_grid(mg, "drainage_area")
figure("K (core only)")
mg.at_node["flow__potential"][mg.boundary_nodes] = 0.
imshow_grid(mg, "flow__potential")
depths = np.where(
    mg.at_node["surface_water__depth"] > 1.e6, 0., mg.at_node["surface_water__depth"]
)
figure("depths")
imshow_grid(mg, depths, var_name="depths (Manning)")
figure("xcomp")
imshow_grid(mg, "water__discharge_x_component")
figure("ycomp")
imshow_grid(mg, "water__discharge_y_component")

print("flux in per node: ", mg.at_node["water__unit_flux_in"][0])
print(
    "water in, water out: ",
    np.sum(mg.at_node["water__unit_flux_in"][mg.core_nodes]),
    np.sum(mg.at_node["surface_water__discharge"][mg.boundary_nodes]),
)
print(
    -np.sum(mg.at_node["surface_water__discharge"][mg.boundary_nodes])
    + np.sum(mg.at_node["water__unit_flux_in"][mg.core_nodes])
)
print(
    (
        -np.sum(mg.at_node["surface_water__discharge"][mg.boundary_nodes])
        + np.sum(mg.at_node["water__unit_flux_in"][mg.core_nodes])
    )
    / mg.at_node["water__unit_flux_in"][0]
)

max_stars_repo_name: cctrunz/landlab

max_issues_repo_name: cctrunz/landlab

max_forks_repo_name: cctrunz/landlab


text:
import numpy as np

from env.sawyer.sawyer import SawyerEnv

class SawyerAssemblyObstacleEnv(SawyerEnv):
    def __init__(self, **kwargs):
        super().__init__("sawyer_assembly_obstacle.xml", **kwargs)
        self._get_reference()

    def _get_reference(self):
        super()._get_reference()

    @property
    def dof(self):
        return 7

    @property
    def init_qpos(self):
        return np.array([0.427, 0.13, 0.0557, 0.114, -0.0622, 0.0276, 0.00356])

    def _reset(self):
        init_qpos = (
            self.init_qpos + self.np_random.randn(self.init_qpos.shape[0]) * 0.02
        )
        self.sim.data.qpos[self.ref_joint_pos_indexes] = init_qpos
        self.sim.data.qvel[self.ref_joint_vel_indexes] = 0.0
        self.sim.forward()

        return self._get_obs()

    def compute_reward(self, action):
        info = {}
        reward = 0
        reward_type = self._kwargs["reward_type"]
        pegHeadPos = self.sim.data.get_site_xpos("pegHead")
        hole = self.sim.data.get_site_xpos("hole")
        dist = np.linalg.norm(pegHeadPos - hole)
        hole_bottom = self.sim.data.get_site_xpos("hole_bottom")
        dist_to_hole_bottom = np.linalg.norm(pegHeadPos - hole_bottom)
        dist_to_hole = np.linalg.norm(pegHeadPos - hole)
        reward_reach = 0
        if dist < 0.3:
            reward_reach += 0.4 * (1 - np.tanh(15 * dist_to_hole))
        reward += reward_reach
        if dist_to_hole_bottom < 0.025:
            reward += self._kwargs["success_reward"]
            self._success = True
            self._terminal = True

        return reward, info

    def _get_obs(self):
        di = super()._get_obs()
        di["hole"] = self.sim.data.get_site_xpos("hole")
        di["pegHead"] = self.sim.data.get_site_xpos("pegHead")
        di["pegEnd"] = self.sim.data.get_site_xpos("pegEnd")
        di["peg_quat"] = self._get_quat("peg")
        return di

    @property
    def static_bodies(self):
        return ["table"]

    @property
    def manipulation_bodies(self):
        return ["furniture", "0_part0", "1_part1", "4_part4", "2_part2"]

    @property
    def manipulation_geom_ids(self):
        body_ids = []
        for body_name in self.manipulation_bodies:
            body_ids.append(self.sim.model.body_name2id(body_name))

        geom_ids = []
        for geom_id, body_id in enumerate(self.sim.model.geom_bodyid):
            if body_id in body_ids:
                geom_ids.append(geom_id)
        return geom_ids

    @property
    def static_geoms(self):
        return []

    @property
    def static_geom_ids(self):
        body_ids = []
        for body_name in self.static_bodies:
            body_ids.append(self.sim.model.body_name2id(body_name))

        geom_ids = []
        for geom_id, body_id in enumerate(self.sim.model.geom_bodyid):
            if body_id in body_ids:
                geom_ids.append(geom_id)
        return geom_ids

    def _step(self, action, is_planner=False):
        """
        (Optional) does gripper visualization after actions.
        """
        assert len(action) == self.dof, "environment got invalid action dimension"

        if not is_planner or self._prev_state is None:
            self._prev_state = self.sim.data.qpos[self.ref_joint_pos_indexes].copy()

        if self._i_term is None:
            self._i_term = np.zeros_like(self.mujoco_robot.dof)

        if is_planner:
            rescaled_ac = np.clip(
                action[: self.robot_dof], -self._ac_scale, self._ac_scale
            )
        else:
            rescaled_ac = np.clip(
                action[: self.robot_dof] * self._ac_scale,
                -self._ac_scale,
                self._ac_scale,
            )
        desired_state = self._prev_state + rescaled_ac

        n_inner_loop = int(self._frame_dt / self.dt)
        for _ in range(n_inner_loop):
            self.sim.data.qfrc_applied[
                self.ref_joint_vel_indexes
            ] = self.sim.data.qfrc_bias[self.ref_joint_vel_indexes].copy()

            if self.use_robot_indicator:
                self.sim.data.qfrc_applied[
                    self.ref_indicator_joint_pos_indexes
                ] = self.sim.data.qfrc_bias[self.ref_indicator_joint_pos_indexes].copy()

            if self.use_target_robot_indicator:
                self.sim.data.qfrc_applied[
                    self.ref_target_indicator_joint_pos_indexes
                ] = self.sim.data.qfrc_bias[
                    self.ref_target_indicator_joint_pos_indexes
                ].copy()
            self._do_simulation(desired_state)

        self._prev_state = np.copy(desired_state)
        reward, info = self.compute_reward(action)

        return self._get_obs(), reward, self._terminal, info

max_stars_repo_name: clvrai/mopa-rl

max_issues_repo_name: clvrai/mopa-rl

max_forks_repo_name: clvrai/mopa-rl


text:
""" Utility functions for returning slice times from number of slices and TR

Slice timing routines in nipy need a vector of slice times.

Slice times are vectors $t_i$ with $i = 0 ... N$ of times, one for each slice, where
$t_i$ gives the time at which slice number $i$ was acquired, relative to the
beginning of the volume acquisition.

We like these vectors because they are unambiguous; the indices $i$ refer to
positions in space, and the values $t_i$ refer to times.

But, there are many common slice timing regimes for which it's easy to get the
slice times once you know the volume acquisition time (the TR) and the number of
slices.

For example, if you acquired the slices in a simple ascending order, and you
have 10 slices and the TR was 2.0, then the slice times are:

>>> import numpy as np
>>> np.arange(10) / 10.  * 2.0
array([ 0. ,  0.2,  0.4,  0.6,  0.8,  1. ,  1.2,  1.4,  1.6,  1.8])

These are small convenience functions that accept the number of slices and the
TR as input, and return a vector of slice times:

>>> ascending(10, 2.)
array([ 0. ,  0.2,  0.4,  0.6,  0.8,  1. ,  1.2,  1.4,  1.6,  1.8])
"""
from __future__ import division, print_function, absolute_import

import numpy as np

# Dictionary (key, value) == (name, function) for slice timing functions
SLICETIME_FUNCTIONS = {}

def _dec_filldoc(func):
    """ Fill docstring of slice time function
    """
    func._doc_template = func.__doc__
    func.__doc__ = func.__doc__.format(
        **dict(
            name = func.__name__,
            pstr=
"""Note: slice 0 is the first slice in the voxel data block

    Parameters
    ----------
    n_slices : int
        Number of slices in volume
    TR : float
        Time to acquire one full volume

    Returns
    -------
    slice_times : (n_slices,) ndarray
        Vectors $t_i i = 0 ... N$ of times, one for each slice, where $t_i$
        gives the time at which slice number $i$ was acquired, relative to the
        beginning of the volume acquisition.
    """))
    return func

def _dec_register_stf(func):
    """ Register slice time function in module dictionary """
    name = func.__name__
    SLICETIME_FUNCTIONS[name] = func
    if name.startswith('st_'):
        short_name = name[3:]
        if short_name in SLICETIME_FUNCTIONS:
            raise ValueError(
                "Duplicate short / long function name {0}".format(short_name))
        SLICETIME_FUNCTIONS[short_name] = func
    return func

def _dec_stfunc(func):
    return _dec_register_stf(_dec_filldoc(func))

def _derived_func(name, func):
    def derived(n_slices, TR):
        return func(n_slices, TR)
    derived.__name__ = name
    derived.__doc__ = func._doc_template
    return _dec_stfunc(derived)

@_dec_stfunc
def st_01234(n_slices, TR):
    """ Simple ascending slice sequence

    slice 0 first, slice 1 second etc.

    For example, for 5 slices and a TR of 1:

    >>> {name}(5, 1.)
    array([ 0. ,  0.2,  0.4,  0.6,  0.8])

    {pstr}
    """
    return np.arange(n_slices) / n_slices * TR

ascending = _derived_func('ascending', st_01234)

@_dec_stfunc
def st_43210(n_slices, TR):
    """ Simple descending slice sequence

    slice ``n_slices-1`` first, slice ``n_slices - 2`` second etc.

    For example, for 5 slices and a TR of 1:

    >>> {name}(5, 1.)
    array([ 0.8,  0.6,  0.4,  0.2,  0. ])

    {pstr}
    """
    return np.arange(n_slices)[::-1] / n_slices * TR

descending = _derived_func('descending', st_43210)

@_dec_stfunc
def st_02413(n_slices, TR):
    """Ascend alternate every second slice, starting at first slice

    Collect slice 0 first, slice 2 second up to top.  Then return to collect
    slice 1, slice 3 etc.

    For example, for 5 slices and a TR of 1:

    >>> {name}(5, 1.)
    array([ 0. ,  0.6,  0.2,  0.8,  0.4])

    {pstr}
    """
    one_slice = TR / n_slices
    time_to_space = list(range(0, n_slices, 2)) + list(range(1, n_slices, 2))
    space_to_time = np.argsort(time_to_space)
    return space_to_time * one_slice

asc_alt_2 = _derived_func('asc_alt_2', st_02413)

@_dec_stfunc
def st_13024(n_slices, TR):
    """Ascend alternate every second slice, starting at second slice

    Collect slice 1 first, slice 3 second up to top (highest numbered slice).
    Then return to collect slice 0, slice 2 etc.  This order is rare except on
    Siemens acquisitions with an even number of slices.  See
    :func:`st_odd0_even1` for this logic.

    For example, for 5 slices and a TR of 1:

    >>> {name}(5, 1.)
    array([ 0.4,  0. ,  0.6,  0.2,  0.8])

    {pstr}
    """
    one_slice = TR / n_slices
    time_to_space = list(range(1, n_slices, 2)) + list(range(0, n_slices, 2))
    space_to_time = np.argsort(time_to_space)
    return space_to_time * one_slice

asc_alt_2_1 = _derived_func('asc_alt_2_1', st_13024)

@_dec_stfunc
def st_42031(n_slices, TR):
    """Descend alternate every second slice, starting at last slice

    Collect slice (`n_slices` - 1) first, slice (`nslices` - 3) second down to
    bottom (lowest numbered slice).  Then return to collect slice (`n_slices`
    -2), slice (`n_slices` - 4) etc.

    For example, for 5 slices and a TR of 1:

    >>> {name}(5, 1.)
    array([ 0.4,  0.8,  0.2,  0.6,  0. ])

    {pstr}
    """
    return st_02413(n_slices, TR)[::-1]

desc_alt_2 = _derived_func('desc_alt_2', st_42031)

@_dec_stfunc
def st_odd0_even1(n_slices, TR):
    """Ascend alternate starting at slice 0 for odd, slice 1 for even `n_slices`

    Acquisitions with alternating ascending slices from Siemens scanners often
    seem to have this behavior as default - see:
        https://mri.radiology.uiowa.edu/fmri_images.html

    This means we use the :func:`st_02413` algorithm if `n_slices` is odd, and
    the :func:`st_13024` algorithm if `n_slices` is even.

    For example, for 4 slices and a TR of 1:

    >>> {name}(4, 1.)
    array([ 0.5 ,  0.  ,  0.75,  0.25])

    5 slices and a TR of 1:

    >>> {name}(5, 1.)
    array([ 0. ,  0.6,  0.2,  0.8,  0.4])

    {pstr}
    """
    if n_slices % 2 == 0:
        return st_13024(n_slices, TR)
    return st_02413(n_slices, TR)

asc_alt_siemens = _derived_func('asc_alt_siemens', st_odd0_even1)

@_dec_stfunc
def st_03142(n_slices, TR):
    """Ascend alternate, where alternation is by half the volume

    Collect slice 0 then slice ``ceil(n_slices / 2.)`` then slice 1 then slice
    ``ceil(nslices / 2.) + 1`` etc.

    For example, for 5 slices and a TR of 1:

    >>> {name}(5, 1.)
    array([ 0. ,  0.4,  0.8,  0.2,  0.6])

    {pstr}
    """
    one_slice = TR / n_slices
    space_to_time = (list(range(0, n_slices, 2)) +
                     list(range(1, n_slices, 2)))
    return np.array(space_to_time) * one_slice

asc_alt_half = _derived_func('asc_alt_half', st_03142)

@_dec_stfunc
def st_41302(n_slices, TR):
    """Descend alternate, where alternation is by half the volume

    Collect slice ``(n_slices - 1)`` then slice ``floor(nslices / 2.) - 1``
    then slice ``(n_slices - 2)`` then slice ``floor(nslices / 2.) - 2`` etc.

    For example, for 5 slices and a TR of 1:

    >>> {name}(5, 1.)
    array([ 0.6,  0.2,  0.8,  0.4,  0. ])

    {pstr}
    """
    return st_03142(n_slices, TR)[::-1]

desc_alt_half = _derived_func('desc_alt_half', st_41302)

max_stars_repo_name: bpinsard/nipy

max_issues_repo_name: bpinsard/nipy

max_forks_repo_name: bpinsard/nipy


text:
import h5py
import sys
sys.path.append("/home/lur159/git/miniMDS")
from tools import Tracker
import numpy as np

in_path = sys.argv[1]
out_path = sys.argv[2]
f = h5py.File(in_path)

counts = np.array(f["pixels"]["count"])
bin_ids1 = np.array(f["pixels"]["bin1_id"])
bin_ids2 = np.array(f["pixels"]["bin2_id"])
chroms = np.array(f["bins"]["chrom"])
starts = np.array(f["bins"]["start"])
ends = np.array(f["bins"]["end"])

f.close()

tracker = Tracker("Converting to BED", len(counts))

print "Begin converting to BED"
with open(out_path, "w") as out_file:
	for count, bin_id1, bin_id2 in zip(counts, bin_ids1, bin_ids2):
		if count != 0:
			chrom1 = str(chroms[bin_id1] + 1)	#switch to 1-indexed
			chrom2 = str(chroms[bin_id2] + 1)
			start1 = str(starts[bin_id1])
			end1 = str(ends[bin_id1])
			start2 = str(starts[bin_id2])
			end2 = str(ends[bin_id2])
			out_file.write("\t".join(("chr" + chrom1, start1, end1, "chr" + chrom2, start2, end2, str(count))))
			out_file.write("\n")
		tracker.increment()
	out_file.close()


text:
import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sp
import scipy.sparse.linalg as la
from functools import partial
import time

# Creating the MAtrix as describve din the report
def make_L(Nx, Ny):
    Dx = sp.diags((Nx-1)*[1])
    Dx += sp.diags((Nx-2)*[-1],-1)
    rowx = sp.csr_matrix((1,Nx-1))
    rowx[0,-1] = -1
    Dx = sp.vstack((Dx, rowx))
    Lx = Dx.transpose().dot(Dx)

    Dy = sp.diags((Ny-1)*[1])
    Dy += sp.diags((Ny-2)*[-1],-1)
    rowy = sp.csr_matrix((1,Ny-1))
    rowy[0,-1] = -1
    Dy = sp.vstack((Dy, rowy))
    Ly = Dy.transpose().dot(Dy)

    #Kronsum is jsut a cleaner way than creating Identity matrices and stuff
    return sp.kronsum(Lx,Ly)

def discretize(x_d, y_d, h):
    nx = int(x_d/h)
    ny = int(y_d/h)
    return make_L(nx,ny)/h/h

def get_grid(x_d, y_d, h):
    grid = np.mgrid[h:y_d:h, h:x_d:h]
    return (grid[1,:,:], grid[0,:,:])

def source(xx,yy):
    a = -10
    b = 5
    return np.exp(a*(np.square(xx-b)+np.square(yy-b)))

def sourcevec(xx,yy):
    return np.reshape(source(xx,yy), (xx.shape[0]*xx.shape[1]))

imshow = partial(plt.imshow)

#Domain is (0,0)x(x,y)
x = 10
y = 10

#choose grid spacing
h = 0.1

gamma = [-40, 0, 40]
solutions = []
sources = []
residues = []

for i in range(3):
    # Creating grid, L
    grid = get_grid(x,y,h)
    L = discretize(x,y,h)
    L = L - gamma[i] * sp.eye(L.shape[0])

    #Creation of the source vector
    sv = sourcevec(*grid)

    #Solving the system
    residuals = []

    def cb(rk):
        print("\rIn iteration number %4d, rk is %1.5e"%(len(residuals)+1,rk), end="")
        residuals.append(rk)
    start = time.time()
    solution, succ = la.gmres(L, sv, maxiter=5000, restart=5000, tol=1e-12, callback=cb)
    # check if GMRES COnverged
    if succ == 0:
        print("\nGMRES Converged")
    elif succ > 0:
        print("\nGMRES Converged but given tolerance not achieved or max iterations reached")
    else:
        print("\nYeah, you made an oopsie")

    print("GMRES took %3.2fs"%(time.time()-start))
    residues.append(residuals)
    solutions.append(solution)
    sources.append(sv)

    print("This should be small:", np.linalg.norm(sv-L@solution)/np.linalg.norm(sv)-residuals[-1])

    start = time.time()
    solution = la.spsolve(L, sv)
    print("spsolve took %3.2fs"%(time.time()-start))

#Showing source then Solution
reshaper = lambda u: np.reshape(u, [grid[0].shape[0], grid[0].shape[1]])[::-1,:]
for i in range(3):
    plt.semilogy(residues[i], label="$\gamma=$%d"%(gamma[i]))
plt.legend()
plt.ylabel("Residual")
plt.xlabel("Iteration")
plt.show()

fig = plt.figure()
plt.subplot(2,2,1)
plt.title("Source function")
plt.imshow(reshaper(sources[0]))
plt.colorbar()
for i in range(3):
    plt.subplot(2,2,i+2)
    plt.title("$\gamma=$%d"%(gamma[i]))
    plt.imshow(reshaper(solutions[i]))
    plt.colorbar()

plt.show()

reason: import numpy,import scipy


text:
from keras import backend as K
from keras.models import Model, model_from_json
from keras.layers import Input, TimeDistributed, Dense, Flatten
import numpy as np
import h5py

from prednet import PredNet

def load_model(model_json_file, model_weights_file, **extras):
    print('Loading model: {}'.format(model_weights_file))
    with open(model_json_file, 'r') as f:
        json_string = f.read()
        train_model = model_from_json(json_string, custom_objects={'PredNet': PredNet})
        train_model.load_weights(model_weights_file)
    return train_model

def create_model(model_json_file=None, model_weights_file=None,
                 train=False, **config):
    if model_json_file and model_weights_file:
        pretrained_model = load_model(model_json_file, model_weights_file)
        model = pretrained_prednet(pretrained_model, train=train, **config)
    else:
        model = random_prednet(train=train, **config)
    return model

def pretrained_prednet(pretrained_model, n_timesteps, output_mode='error',
                       train=False, stateful=False, batch_size=None,
                       trainable_layers=None, trainable_units=None,
                       **config):
    prednet_model = pretrained_model

    if 'stack_sizes' not in prednet_model.layers[1].get_config():
        for layer in pretrained_model.layers:
            if 'prednet' in layer.name.lower():
                print('Found PredNet in layer', layer.name)
                prednet_model = layer
                break

    layer_config = prednet_model.layers[1].get_config()
    layer_config['output_mode'] = output_mode
    layer_config['stateful'] = stateful
    prednet = PredNet(weights=prednet_model.layers[1].get_weights(),
                      trainable_layers=trainable_layers,
                      trainable_units=trainable_units, **layer_config)
    input_shape = list(prednet_model.layers[0].batch_input_shape[1:])
    input_shape[0] = n_timesteps
    inputs = get_input_layer(batch_size, tuple(input_shape))
    outputs = get_output_layer(prednet, inputs, n_timesteps, train, output_mode)
    model = Model(inputs=inputs, outputs=outputs, name='PredNet')
    return model

def random_prednet(input_channels, input_height, input_width,
                   n_timesteps, stack_sizes=(48, 96, 192),
                   train=False, output_mode='error', stateful=False,
                   batch_size=None, trainable_layers=None,
                   trainable_units=None, **config):
    # Model parameters
    if K.image_data_format() == 'channels_first':
        input_shape = (input_channels, input_height, input_width)
    else:
        input_shape = (input_height, input_width, input_channels)

    stack_sizes = (input_channels,) + stack_sizes
    R_stack_sizes = stack_sizes
    A_filt_sizes = (3,) * (len(stack_sizes) - 1)
    Ahat_filt_sizes = (3,) * len(stack_sizes)
    R_filt_sizes = (3,) * len(stack_sizes)
    prednet = PredNet(stack_sizes, R_stack_sizes,
                      A_filt_sizes, Ahat_filt_sizes, R_filt_sizes,
                      output_mode=output_mode, return_sequences=True,
                      stateful=stateful, trainable_layers=trainable_layers,
                      trainable_units=trainable_units, name='prednet_layer')
    input_shape = (n_timesteps,) + input_shape
    inputs = get_input_layer(batch_size, input_shape)
    outputs = get_output_layer(prednet, inputs, n_timesteps, train, output_mode)
    model = Model(inputs=inputs, outputs=outputs, name='PredNet')
    return model

def get_input_layer(batch_size, input_shape):
    if batch_size:
        input_shape = (batch_size,) + input_shape
        inputs = Input(batch_shape=input_shape)
    else:
        inputs = Input(shape=input_shape)
    return inputs

def get_output_layer(prednet, inputs, n_timesteps, train, output_mode):
    outputs = prednet(inputs)
    if train:
        if output_mode not in ['error', 'representation_and_error']:
            raise ValueError('When training, output_mode must be equal to "error"')
    if output_mode == 'error':
        outputs = get_error_layer(outputs, n_timesteps, prednet.nb_layers)
    return outputs

def get_error_layer(outputs, n_timesteps, nb_layers):
    # weighting for each layer in final loss; "L_0" model:  [1, 0, 0, 0], "L_all": [1, 0.1, 0.1, 0.1]
    layer_loss_weights = np.array([0.1] * nb_layers)
    layer_loss_weights[0] = 1
    layer_loss_weights = np.expand_dims(layer_loss_weights, 1)
    # equally weight all timesteps except the first
    time_loss_weights = 1. / (n_timesteps - 1) * np.ones((n_timesteps, 1))
    time_loss_weights[0] = 0

    # calculate weighted error by layer
    errors_by_time = TimeDistributed(Dense(1, trainable=False),
                                     weights=[layer_loss_weights, np.zeros(1)],
                                     trainable=False)(outputs)
    errors_by_time = Flatten()(errors_by_time)  # will be (batch_size, n_timesteps)
    final_errors = Dense(1, weights=[time_loss_weights, np.zeros(1)],
                         trainable=False)(errors_by_time)  # weight errors by time
    return final_errors

def load_prednet_weights(weights_path, prednet_model):
    """
    Load model weights from HDF5 directly and set weights tensor by tensor.
    This implementation is specific for a PredNet model.
    This is workaround for some issues related to HDF5 loading on Google Colab:
    OSError: Unable to open file (File signature not found)

    :param prednet_model: Keras model instance to set the weights
    :param weights_path: path for HDF5 weights file
    """

    weights = h5py.File(weights_path)
    weight_keys = list(weights.keys())
    if 'model_weights' in weight_keys:
        weights = weights['model_weights']
        weight_keys = list(weights.keys())

    prednet_key = None
    for key in weight_keys:
        if 'prednet' in key.replace('_', ''):
            prednet_key = key
            break

    if prednet_key is None:
        raise ValueError('Weights for prednet layer not found')

    prednet_weights_keys = list(weights[prednet_key].keys())
    prednet_weights = weights[prednet_key][prednet_weights_keys[0]]

    total_weights = 0
    weight_value_tuples = []
    for layer_key in prednet_weights:
        layer_weights = prednet_weights[layer_key]
        total_weights += len(layer_weights)

        for weight_key in layer_weights:
            # print(f'Loading weights for layer {layer_key}, {weight_key} shape = {layer_weights[weight_key].shape}')
            tensor_name = f'{layer_key}/{weight_key}'

            for tensor in prednet_model.layers[1].weights:
                if tensor_name in tensor.name:
                    # print(f'Setting tensor {tensor.name} to value in weights[{tensor_name}] - shape = {layer_weights[weight_key].shape}')
                    weight_value_tuples.append((tensor, layer_weights[weight_key]))

    if weight_value_tuples:
        K.batch_set_value(weight_value_tuples)
    print(f'{len(weight_value_tuples)}/{total_weights} weights loaded from {weights_path}')

def load_model_weights(weights_path, target_layer, source_weights_key=None):
    """
    Load model weights from HDF5 directly and set weights tensor by tensor.
    This is workaround for some issues related to HDF5 loading on Google Colab:
    OSError: Unable to open file (File signature not found)

    :param weights_path: path for HDF5 weights file
    :param target_layer: layer for which we want to set the weights
    :param source_weights_key: some specific HDF5 key for the weights
    """
    weights = h5py.File(weights_path)
    weight_keys = list(weights.keys())
    if 'model_weights' in weight_keys:
        weights = weights['model_weights']
        weight_keys = list(weights.keys())

    if source_weights_key is None:
        source_weights_key = target_layer.name

    model_weights_key = None
    for key in weight_keys:
        if source_weights_key in key:
            model_weights_key = key
            break

    if model_weights_key is None:
        raise ValueError(f'Weights for layer {target_layer.name} not found')

    layer_weights_keys = list(weights[model_weights_key].keys())
    all_layer_weights = weights[model_weights_key][layer_weights_keys[0]]

    total_weights = 0
    weight_value_tuples = []
    for layer_key in all_layer_weights:
        layer_weights = all_layer_weights[layer_key]
        total_weights += len(layer_weights)

        if type(layer_weights[0]) == str:
            for weights in layer_weights:
                # print(f'Loading weights for layer {layer_key}, {weights} shape = {layer_weights[weights].shape}')
                tensor_name = f'{layer_key}/{weights}'
                weights = layer_weights[weights]

                for tensor in target_layer.weights:
                    print(tensor_name, tensor.name)
                    if tensor_name in tensor.name:
                        weight_value_tuples.append((tensor, weights))
                        break

        else:
            weights = layer_weights
            # print(f'Loading weights for layer {layer_key} shape = {weights.shape}')
            tensor_name = f'{layer_key}'

            for tensor in target_layer.weights:
                print(tensor_name, tensor.name)
                if tensor_name in tensor.name:
                    # print(f'Setting tensor {tensor.name} to value in weights[{tensor_name}] - shape = {weights.shape}')
                    weight_value_tuples.append((tensor, weights))
                    break

    if weight_value_tuples:
        K.batch_set_value(weight_value_tuples)
    print(f'{len(weight_value_tuples)}/{total_weights} weights loaded from {weights_path}')

max_stars_repo_name: thefonseca/algonauts

max_issues_repo_name: thefonseca/algonauts

max_forks_repo_name: thefonseca/algonauts