text:
import numpy as np
import pytest
from sklearn.model_selection import cross_validate

from tslearn.utils import to_time_series, to_time_series_dataset

__author__ = 'Romain Tavenard romain.tavenard[at]univ-rennes2.fr'

def test_shapelets():
    pytest.importorskip('tensorflow')
    from tslearn.shapelets import LearningShapelets
    import tensorflow as tf

    n, sz, d = 15, 10, 2
    rng = np.random.RandomState(0)
    time_series = rng.randn(n, sz, d)
    y = rng.randint(2, size=n)
    clf = LearningShapelets(n_shapelets_per_size={2: 5},
                            max_iter=1,
                            verbose=0,
                            optimizer="sgd",
                            random_state=0)

    cross_validate(clf, time_series, y, cv=2)

    clf = LearningShapelets(n_shapelets_per_size={2: 5},
                            max_iter=1,
                            verbose=0,
                            optimizer=tf.optimizers.Adam(.1),
                            random_state=0)
    cross_validate(clf, time_series, y, cv=2)

    model = LearningShapelets(n_shapelets_per_size={3: 2, 4: 1}, max_iter=1)
    model.fit(time_series, y)
    for shp, shp_bis in zip(model.shapelets_,
                            model.shapelets_as_time_series_):
        np.testing.assert_allclose(shp,
                                   to_time_series(shp_bis, remove_nans=True))

    # Test set_weights / get_weights
    clf = LearningShapelets(n_shapelets_per_size={2: 5},
                            max_iter=1,
                            verbose=0,
                            random_state=0)
    clf.fit(time_series, y)
    preds_before = clf.predict_proba(time_series)
    weights = clf.get_weights()
    # Change number of iterations, then refit, then set weights
    clf.max_iter *= 2
    clf.fit(time_series, y)
    clf.set_weights(weights)
    np.testing.assert_allclose(preds_before,
                               clf.predict_proba(time_series))

def test_shapelet_lengths():
    pytest.importorskip('tensorflow')
    from tslearn.shapelets import LearningShapelets

    # Test variable-length
    y = [0, 1]
    time_series = to_time_series_dataset([[1, 2, 3, 4, 5], [3, 2, 1]])
    clf = LearningShapelets(n_shapelets_per_size={3: 1},
                            max_iter=1,
                            verbose=0,
                            random_state=0)
    clf.fit(time_series, y)

    weights_shapelet = [np.array([[1, 2, 3]])]
    clf.set_weights(weights_shapelet, layer_name="shapelets_0_0")
    tr = clf.transform(time_series)
    np.testing.assert_allclose(tr,
                               np.array([[0.], [8. / 3]]))

    # Test max_size to predict longer series than those passed at fit time
    y = [0, 1]
    time_series = to_time_series_dataset([[1, 2, 3, 4, 5], [3, 2, 1]])
    clf = LearningShapelets(n_shapelets_per_size={3: 1},
                            max_iter=1,
                            verbose=0,
                            max_size=6,
                            random_state=0)
    clf.fit(time_series[:, :-1], y)  # Fit with size 4
    weights_shapelet = [np.array([[1, 2, 3]])]
    clf.set_weights(weights_shapelet, layer_name="shapelets_0_0")
    tr = clf.transform(time_series)
    np.testing.assert_allclose(tr,
                               np.array([[0.], [8. / 3]]))

def test_series_lengths():
    pytest.importorskip('tensorflow')
    from tslearn.shapelets import LearningShapelets

    # Test long shapelets
    y = [0, 1]
    time_series = to_time_series_dataset([[1, 2, 3, 4, 5], [3, 2, 1]])
    clf = LearningShapelets(n_shapelets_per_size={8: 1},
                            max_iter=1,
                            verbose=0,
                            random_state=0)
    np.testing.assert_raises(ValueError, clf.fit, time_series, y)

    # Test small max_size
    y = [0, 1]
    time_series = to_time_series_dataset([[1, 2, 3, 4, 5], [3, 2, 1]])
    clf = LearningShapelets(n_shapelets_per_size={3: 1},
                            max_iter=1,
                            verbose=0,
                            max_size=4,
                            random_state=0)
    np.testing.assert_raises(ValueError, clf.fit, time_series, y)

max_stars_repo_name: andriyor/tslearn

max_issues_repo_name: andriyor/tslearn

max_forks_repo_name: andriyor/tslearn


text:
from utility import *
import random
import numpy as np

@maindeco
def main_menue():

    message_display('Hauptmenü', display_width*0.5, display_width*0.15, 58 )

    button( display_width*0.34375, display_height*0.4, display_width*0.3125,display_height/6,'Training', second_menue)
    button( display_width*0.34375, display_height*0.6, display_width*0.3125,display_height/6,'Training', second_menue)

@maindeco
def second_menue():

    button( display_width*0.1, display_height*0.3, display_width*0.3125,display_height/6,'+', add.display_task, size = 60)
    button( display_width*0.55, display_height*0.3, display_width*0.3125,display_height/6,'-',  main_menue, size = 100)
    button( display_width*0.1, display_height*0.5, display_width*0.3125,display_height/6,'*', main_menue, size = 60)
    button( display_width*0.55, display_height*0.5, display_width*0.3125,display_height/6,'/', main_menue, size = 60)

class base_task(object):

    size = 50
    Input = label(display_width*0.7, display_height*0.5, size = 40)

    def __init__(self,operations, task_num = 10, mixed = False ):

        self.task_num = task_num
        self.mixed = mixed
        #a list containing operator dict keys
        self.operations = operations
        self.input = label(display_width*0.7, display_height*0.5, size = 40)
        self.record = []
        self.current_task = 0

        self.init()

    def generate_numbers(self, start = 0, stop = 20):

        self.nums = []

        kardinal = len(self.operations) +1

        for i in range(kardinal):
            self.nums.append(random.randint(start,stop))

        return self.nums

    def operation_to_string(self):

        self.term = ''.join(str(a) + self.operations[i] for i,a in enumerate(self.nums[:-1])) + str(self.nums[-1])
        self.equation =  self.term + ' ='
        return self.equation

    @maindeco
    def display_task(self):

        message_display(self.operation_to_string(),display_width/2,display_height/2, fix_right = 1)
        self.input.fill_values()
        self.input.display()
        button(display_width*0.55, display_height*0.9, display_width*0.3125,display_height/6
               ,'Eingabe', self.enter, size = 50)

    def get_result(self):

            return eval(self.term)

    def verify(self):

       if self.Input.number == self.get_result():

           return 1
       else:
           return 0

    def enter(self):

        self.record.append(self.verify())

        if self.current_task < self.task_num:
            self.Input.delete()
            self.generate_numbers()
            self.current_task += 1
            self.operation_to_string
        else:
            main_menue()
    def init(self):

        self.generate_numbers()
        self.operation_to_string()

add = base_task(['+'])

ext: py

reason: import numpy


text:
from math import sqrt
import numpy as np
from pwtools.crys import coord_trans, coord_trans3d
from pwtools import num
from pwtools.test.tools import aaae
rand = np.random.rand

def test_coord_trans():
    #-----------------------------------------------------------
    # 2D
    #-----------------------------------------------------------
    c_X = rand(20,3)
    # basis vecs are assumed to be rows
    X = rand(3,3)*5
    Y = rand(3,3)*3

    # transform and back-transform
    c_Y = coord_trans(c_X, old=X, new=Y)
    c_X2 = coord_trans(c_Y, old=Y, new=X)
    aaae(c_X, c_X2)

    # simple dot product must produce same cartesian results:
    # X . v_X = I . v_I = v_I
    X = np.identity(3)
    Y = rand(3,3)*3
    c_X = rand(20,3)
    c_Y = coord_trans(c_X, old=X, new=Y)
    # normal back-transform
    c_X2 = coord_trans(c_Y, old=Y, new=X)
    # 2 forms w/ dot(), assume: basis vecs = rows of X and Y
    c_X3 = np.dot(c_Y, Y)
    c_X4 = np.dot(Y.T, c_Y.T).T
    aaae(c_X, c_X2)
    aaae(c_X, c_X3)
    aaae(c_X, c_X4)

    # some textbook example
    #
    v_I = np.array([1.0,1.5])
    I = np.identity(2)
    # basis vecs as rows
    X = sqrt(2)/2.0*np.array([[1,-1],[1,1]]).T
    Y = np.array([[1,1],[0,1]]).T

    # "identity" transform
    aaae(coord_trans(v_I,I,I), v_I)

    # v in basis X and Y
    v_X = coord_trans(v_I,I,X)
    v_Y = coord_trans(v_I,I,Y)
    aaae(v_X, np.array([1.76776695, 0.35355339]))
    aaae(v_Y, np.array([-0.5,  1.5]))

    # back-transform
    aaae(coord_trans(v_X,X,I), v_I)
    aaae(coord_trans(v_Y,Y,I), v_I)

    # higher "x,y,z"-dims: 4-vectors
    c_X = rand(20,4)
    X = rand(4,4)*5
    Y = rand(4,4)*3
    c_Y = coord_trans(c_X, old=X, new=Y)
    c_X2 = coord_trans(c_Y, old=Y, new=X)
    aaae(c_X, c_X2)

    #-----------------------------------------------------------
    # 3D
    #-----------------------------------------------------------

    # x,y,z case
    c_X = rand(20,3,10)
    X = rand(3,3)*5
    Y = rand(3,3)*3
    c_Y = coord_trans(c_X, old=X, new=Y, axis=1)
    c_X2 = coord_trans(c_Y, old=Y, new=X, axis=1)
    aaae(c_X, c_X2)

    c_X = rand(20,10,3)
    c_Y = coord_trans(c_X, old=X, new=Y, axis=-1)
    c_X2 = coord_trans(c_Y, old=Y, new=X, axis=-1)
    aaae(c_X, c_X2)

    c_X = rand(3,20,10)
    c_Y = coord_trans(c_X, old=X, new=Y, axis=0)
    c_X2 = coord_trans(c_Y, old=Y, new=X, axis=0)
    aaae(c_X, c_X2)

    # 3d, higher "x,y,z"-dims, i.e. 4-vectors: trajectory of 5 atoms, 10 steps,
    # "4d-coordinates"
    c_X = rand(20,4,10)
    X = rand(4,4)*5
    Y = rand(4,4)*3
    c_Y = coord_trans(c_X, old=X, new=Y, axis=1)
    c_X2 = coord_trans(c_Y, old=Y, new=X, axis=1)
    aaae(c_X, c_X2)

    #-----------------------------------------------------------
    # ND
    #-----------------------------------------------------------

    # arbitrary collection of 4-vectors
    c_X = rand(20,4,10,8)
    X = rand(4,4)*5
    Y = rand(4,4)*3
    c_Y = coord_trans(c_X, old=X, new=Y, axis=1)
    c_X2 = coord_trans(c_Y, old=Y, new=X, axis=1)
    aaae(c_X, c_X2)

    #-----------------------------------------------------------
    # special case 3d
    #-----------------------------------------------------------
    # Note that axis=1 is always the xyz-axis (length 3) if the timeaxis
    # (length 10) would be removed from all arrays (2d case then).
    c_X = rand(20,3,10)
    X = rand(3,3,10)*5
    Y = rand(3,3,10)*3
    c_Y = coord_trans3d(c_X, old=X, new=Y, axis=1, timeaxis=2)
    c_X2 = coord_trans3d(c_Y, old=Y, new=X, axis=1, timeaxis=2)
    aaae(c_X, c_X2)

    c_X = rand(20,10,3)
    X = rand(3,10,3)*5
    Y = rand(3,10,3)*3
    c_Y = coord_trans3d(c_X, old=X, new=Y, axis=1, timeaxis=1)
    c_X2 = coord_trans3d(c_Y, old=Y, new=X, axis=1, timeaxis=1)
    aaae(c_X, c_X2)

    c_X = rand(10,20,3)
    X = rand(10,3,3)*5
    Y = rand(10,3,3)*3
    c_Y = coord_trans3d(c_X, old=X, new=Y, axis=1, timeaxis=0)
    c_X2 = coord_trans3d(c_Y, old=Y, new=X, axis=1, timeaxis=0)
    aaae(c_X, c_X2)

max_stars_repo_name: elcorto/pwtools

max_issues_repo_name: elcorto/pwtools

max_forks_repo_name: elcorto/pwtools


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

# # (baseline development) Silicon per M2 Calculations

# This journal documents the methods and assumptions made to create a baseline material file for silicon.

# ## Mass per M2

# The mass of silicon contained in a PV module is dependent on the size, thickness and number of cells in an average module. Since there is a range of sizes and number of cells per module, we will attempt a weighted average. These weighted averages are based on ITRPV data, which goes back to 2010, Fraunhofer data back to 1990, and

# In[1]:

import numpy as np
import pandas as pd
import os,sys
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 22})
plt.rcParams['figure.figsize'] = (12, 8)

density_si = 2.3290 #g/cm^3 from Wikipedia of Silicon (https://en.wikipedia.org/wiki/Silicon)
#it might be better to have mono-Si and multi-Si densities, including dopants,
#but that density is not readily available

# A Fraunhofer report indicates that in 1990, wafers were 400 micron thick, decreasing to the more modern 180 micron thickness by 2008. ITRPVs back to 2010 indicate that 156 mm x 156mm was the standard size wafer through 2015.

# In[2]:

#now lets try to do this for 2019 through 2030 all at once with dataframes
#taking the average of the ranges specified in ITRPVs

#first we input the market share data for mcSi and monoSi, read in from csv
cwd = os.getcwd() #grabs current working directory
skipcols=['Source']
mrktshr_cellsize = pd.read_csv(cwd+"/../../PV_ICE/baselines/SupportingMaterial/MarketShare_CellSize.csv",
                               index_col='Year', usecols=lambda x: x not in skipcols)
mrktshr_cellsize /=100 #turn whole numbers into decimal percentages
#print(mrktshr_cellsize)

#then split them into two dataframes for later computations
dfmarketshare_mcSi = mrktshr_cellsize.filter(regex = 'mcSi')
dfmarketshare_monoSi = mrktshr_cellsize.filter(regex = 'monoSi')
#adjust column names for matching computation later
dfmarketshare_mcSi.columns = ['share156','share156.75','share157.75','share163','share166up']
dfmarketshare_monoSi.columns = ['share156','share156.75','share157.75','share163','share166up']

print(dfmarketshare_mcSi)
print(dfmarketshare_monoSi)

# Interpolate marketshare for missing years in ITRPV 2020 predictions
# ----
# choosing to interpolate market share of different sizes rather than cell size because this should be more basedin technology - i.e. crystals only grow certain sizes. Additionally, it is more helpful to understand the impact silicon usage by keeping cell size and marketshare seperate.

# In[3]:

#interpolate for missing marketshare data
##the interpolate function returns a view of the df, doesn't modify the df itself
##therefore you have to set the old df, or a new one = the df.interpolate function
dfmarketshare_mcSi=dfmarketshare_mcSi.interpolate(method='linear',axis=0,limit=2,limit_area='inside')
dfmarketshare_monoSi=dfmarketshare_monoSi.interpolate(method='linear',axis=0,limit=2,limit_area='inside')

#fill remaining NaN/outside with 0 (i.e., no market share)
dfmarketshare_mcSi=dfmarketshare_mcSi.fillna(0.0)
dfmarketshare_monoSi=dfmarketshare_monoSi.fillna(0.0)

print(dfmarketshare_mcSi)
print(dfmarketshare_monoSi)

# In[4]:

#multiply each marketshare dataframe column by it's respective size
#dfmarketshare_mcSi.share156 *=156 #this is a slow way to multiply each column by its respective size

cellsizes = {'share156':156,
            'share156.75':156.75,
            'share157.75':157.75,
            'share163':163.875,
            'share166up':166} #dictionary of the average cell dimension for each market share bin (ITRPV 2020)

#multiply cell dimensions by their market share to get a weighted average
##this is where the column names needed to match
df_scalecell_mcSi = dfmarketshare_mcSi.mul(cellsizes,'columns')
df_scalecell_monoSi = dfmarketshare_monoSi.mul(cellsizes,'columns')

print(df_scalecell_mcSi)
print(df_scalecell_monoSi)

# In[5]:

#now add the columns together to get the weighted average cell size for each year for each technology
df_avgcell_mcSi = pd.DataFrame(df_scalecell_mcSi.agg("sum", axis="columns"))
df_avgcell_monoSi = pd.DataFrame(df_scalecell_monoSi.agg("sum", axis="columns")) #agg functions return a series not a dictionary
#print(df_avgcell_mcSi)

#join the two dataframes into single one with two columns
df_avgcell = pd.concat([df_avgcell_monoSi,df_avgcell_mcSi], axis=1) #concatinate on the columns axis
df_avgcell.columns = ['monoSi','mcSi'] #name the columns
print(df_avgcell)

# However, we know that it wasn't 156 mm back to 1995, but exact records of average cell size are lacking. A mention of a companies' new manufacturing line producing 125 mm mono-Si in 1993 can be found in IEA PVPS documentation, and Martin Green 2000 calls out 100 mm to 150 mm manufacturing. SinoVoltaics "History of Wafer Sizes" claims 100 mm x 100mm was standard through 1996, then 125 mm x 125mm took over, followed by 156 mm (start date not specified), with 156.75mm starting in 2016. A 2020 Trina solar presentation claims 100 mm x 100 mm in 1995, and only increases to 125mm in 2005, with 156mm taking over in ~2010 (which generally agrees with ITRPVs).
#
# WoodMackenzie Reports indicate that 1981 through 2012, 125mm was the dominant size, 2013-2016 was 156mm and 156.75mm, 2017-2018 had many sizes between 158.75mm and 161.75mm, and finally 2019 saw the introduction of 166m, 182mm and 210mm (Sun, Solar PV Module Technology Report 2020). The report data notes that these are "shipment" or manufacturing data, and the installed capacity will lag by 3 to 18 months. Finally, the report predicts that 182mm or 210mm will become dominant shipped technology by 2023. This timeline is ahead of what the ITRPVs predict, though may align when accounting for the 3 to 18 month lag.
#
# Based on these sources, we will say that cell sizes in 1995 were 100 mm, and in 2000 were 125 mm, and 156mm in 2010 (where ITRPV data starts). These will be step functions instead of linear interpolations to attempt to better represent that most size changes requires a replacement of the manufacturing line equipment. In reality, there would be some marketshare blending, which could be added with improved future data finding.

# In[6]:

#turn zeros back into NaN
df_avgcell.replace(0.0, np.NaN, inplace=True)

#write over 1995 and 2000 data
#df_avgcell['monoSi'][1995]=100.00
#df_avgcell['monoSi'][2000]=125.00
#df_avgcell['mcSi'][1995]=100.00
#df_avgcell['mcSi'][2000]=125.00

#set dates for start and end of cell size changes
start_100 = 1995
end_100 = 1999
start_125 = 2000
end_125 = 2009

df_avgcell.loc[start_100:end_100] = 100.00
df_avgcell.loc[start_125:end_125] = 125.00

print(df_avgcell)

# Now we have an average cell dimension for mc-Si and mono-Si for 1995 through 2030.

# ## Marketshare Data Manipulation

# Next, we apply the marketshare of mc-Si vs mono-Si to get the average cell dimension for the year. Market share of mc-Si vs mono-Si is taken from LBNL "Tracking the Sun" report (warning: this is non-utility scale data i.e. <5MW, and is from 2002-2018), from Mints 2019 SPV report, from ITRPVs, and old papers (Costello & Rappaport 1980, Maycock 2003 & 2005).

# In[7]:

#read in a csv that was copied from CE Data google sheet
cwd = os.getcwd() #grabs current working directory
techmarketshare = pd.read_csv(cwd+"/../../PV_ICE/baselines/SupportingMaterial/ModuleType_MarketShare.csv",index_col='Year')
#this file path navigates from current working directory back up 2 folders, and over to the csv
techmarketshare /=100 #turn whole numbers into decimal percentages
print(techmarketshare)

# #### create a harmonization of annual market share, and interpolate

# In[8]:

# first, create a single value of tech market share in each year or NaN
#split mcSi and monoSi
mcSi_cols = techmarketshare.filter(regex = 'mcSi')
monoSi_cols = techmarketshare.filter(regex = 'mono')

#show where the data is coming from graphically
#plt.plot(mcSi_cols,'.')
monoSikeys = monoSi_cols.columns
labelnames = [e[7:] for e in monoSikeys] #e is a random variable to iterate through the elements of the array
mcSikeys = mcSi_cols.columns
labelnames_mcSi = [e[5:] for e in mcSikeys]
#print(monoSikeys)

# In[9]:

#aggregate all the columns of mono or mcSi into one averaged market share
est_mktshr_mcSi = pd.DataFrame(mcSi_cols.agg("mean", axis="columns"))
#print(est_mktshr_mcSi)
est_mktshr_monoSi = pd.DataFrame(monoSi_cols.agg("mean", axis="columns"))
#print(est_mktshr_monoSi)

#Join the monoSi and mcSi back together as a dataframe
est_mrktshrs = pd.concat([est_mktshr_monoSi,est_mktshr_mcSi], axis=1) #concatinate on the columns axis
est_mrktshrs.columns = ['monoSi','mcSi'] #name the columns

#plot individuals AND average aggregate
plt.plot(monoSi_cols.index,monoSi_cols[monoSikeys[0]],lw=2,marker='o',label=labelnames[0])
plt.plot(monoSi_cols.index,monoSi_cols[monoSikeys[1]],lw=2,marker='D',label=labelnames[1])
plt.plot(monoSi_cols.index,monoSi_cols[monoSikeys[2]],lw=2,marker='P',label=labelnames[2])
plt.plot(monoSi_cols.index,monoSi_cols[monoSikeys[3]],lw=2,marker='^',label=labelnames[3])
plt.plot(monoSi_cols.index,monoSi_cols[monoSikeys[4]],lw=2,marker='o',label=labelnames[4])
plt.plot(est_mrktshrs.index,est_mrktshrs['monoSi'],'--k', marker='s',label='avgd marketshare of mono')
plt.legend()
plt.xlim([1979,2031])
plt.title('monoSi Market Share by Source')
plt.xlabel('Year')
plt.ylabel('Market Share (%)')

# In[10]:

plt.plot(mcSi_cols.index,mcSi_cols[mcSikeys[0]],lw=2,marker='o',label=labelnames_mcSi[0])
plt.plot(mcSi_cols.index,mcSi_cols[mcSikeys[1]],lw=2,marker='D',label=labelnames_mcSi[1])
plt.plot(mcSi_cols.index,mcSi_cols[mcSikeys[2]],lw=2,marker='P',label=labelnames_mcSi[2])
plt.plot(mcSi_cols.index,mcSi_cols[mcSikeys[3]],lw=2,marker='^',label=labelnames_mcSi[3])
plt.plot(est_mrktshrs.index,est_mrktshrs['mcSi'],'--k', marker='s',label='avgd marketshare of mcSi')
plt.legend()
plt.title('mcSi Market Share by Source')
plt.xlabel('Year')
plt.ylabel('Market Share (%)')

# ### Interpolate and Normalize

# In[11]:

#Interpolate for marketshare NaN values
est_mrktshrs['mcSi'][1980]=0.0
est_mrktshrs = est_mrktshrs.interpolate(method='linear',axis=0,limit_area='inside')

#sanity check of market share data - does it add up?
est_mrktshrs['Total'] = est_mrktshrs.monoSi+est_mrktshrs.mcSi
plt.plot(est_mrktshrs['monoSi'], label='MonoSi')
plt.plot(est_mrktshrs['mcSi'], label='mcSi')
plt.plot(est_mrktshrs['Total'],label='Total')
plt.legend()
plt.title('Annual Technology Market Share')
plt.xlabel('Year')
plt.ylabel('Market Share (%)')
#print(est_mrktshrs)
#Warning: 2002, 10% of the silicon marketshare was "other", including amorphous, etc.
#del est_mrktshrs['Total']

# In[12]:

#normalize all marketshares each year to make sure everything adds to 100%
est_mrktshrs['Scale'] = 1/est_mrktshrs['Total']
est_mrktshrs['monoSi_scaled']= est_mrktshrs['Scale']*est_mrktshrs['monoSi']
est_mrktshrs['mcSi_scaled']= est_mrktshrs['Scale']*est_mrktshrs['mcSi']

scaled_marketshares = est_mrktshrs[['monoSi_scaled','mcSi_scaled']]
scaled_marketshares.columns = ['monoSi','mcSi']
scaled_marketshares.to_csv(cwd+'/../../PV_ICE/baselines/SupportingMaterial/output_scaledmrktshr_mcSi_mono.csv', index=True)
scaled_marketshares['Total'] = scaled_marketshares['monoSi']+scaled_marketshares['mcSi']
#print(scaled_marketshares)
plt.plot(scaled_marketshares['monoSi'],label='MonoSi')
plt.plot(scaled_marketshares['mcSi'],label='mcSi')
plt.plot(scaled_marketshares['Total'],label='Total')
plt.legend()
plt.title('Annual Technology Market Share - normalized')
plt.xlabel('Year')
plt.ylabel('Market Share (%)')

# Combining Cell Size shares and Market Share Data
# ----------
# Now we have separate mono and mcSi dataframes, which contain the average cell size, based on the market share of the cell size bin as enumerated in ITRPV 2020. The next step is to combine these technology specific (mono vs mc) based on the module technology market share.

# In[13]:

#now combine technology market share of mcSi and monoSi with their respective average cell dimensions
#which have already been cell size marketshare weighted
#going to ignore "otherSi" because for the most part less than 2%, except 2002

#trim the techmarketshare data to 1995 through 2030
est_mrktshrs_sub = scaled_marketshares.loc[est_mrktshrs.index>=1995] #could also use a filter function instead

#multiply the share of each tech by the weighted average cell size
mrkt_wtd_cells = est_mrktshrs_sub.mul(df_avgcell,'columns')
#sum across monoSi and mcSi for the total market average cell size (x and y)
market_average_cell_dims = pd.DataFrame(mrkt_wtd_cells.agg("sum", axis="columns"))
market_average_cell_dims.columns = ['avg_cell']

#print(market_average_cell_dims)
plt.plot(market_average_cell_dims, label='annual average cell dimensions in mm')
#plt.legend()
plt.title('Average cell dimensions each year')
plt.xlabel('Year')
plt.ylabel('Average cell dimension (mm)')

# Area of a cell
# -------
# The above weighted averages are 1 axis dimension of the square cells in a module. Here we create a dataframe of the averge area of a cell for each year.

# In[14]:

df_cellarea_mm2 = market_average_cell_dims.pow(2,axis='columns') #still in mm^2/cell
#you cannot multiply the df.columnname by itself and get a dataframe back, but df.pow returns a dataframe
df_cellarea_mm2.columns = ['avg_cell']
df_cellarea_m2 = df_cellarea_mm2*0.000001 #mm^2 to m^2
df_cellarea_cm2 = df_cellarea_mm2*0.01 #mm^2 to cm^2
#print(df_cellarea_m2)

# ## Calculate cells/m^2

# While there is technology information of # of cells per module (ex: 60, 72), we are looking for the amount of silicon per m^2 of module. Therefore, it will be easier to figure out how many cells at their respective sizes fit into a m^2 rather than scaling up to the module, only to divide by module size. Additionally, the analysis excludes any spacing efficiency factor (i.e. how close the cells are together), as this type of information is not readily available. Therefore, the assumption is the cells are placed close together, leaving no space, which should slightly overestimate the silicon per m^2.
#
# This # cells/ m^2 of module will be used as a factor in the final calculation of g Si/m^2 module.

# In[15]:

# calculate # cells/m^2 at this point, rather than using the # cells per module factor
df_cellperm2 = 1/df_cellarea_m2
#print(df_cellperm2)
df_cellperm2.to_csv(cwd+'/../../PV_ICE/baselines/SupportingMaterial/output_cell_per_m2.csv', index=True)

# g of Si per cell
# ---------
# In addition to the number of cells that fit into 1m^2 of module, we need the weight of silicon per cell. First, the weighted average of wafer thickness was calculated for each year based on wafer trends and module type market share in the CE Data google spreadsheet. This data is read in here.

# In[16]:

#read in a csv that was copied from CE Data google sheet where the marketshare weighting was done
cwd = os.getcwd() #grabs current working directory
wafer_thickness = pd.read_csv(cwd+"/../../PV_ICE/baselines/SupportingMaterial/Wafer_thickness.csv",
                              index_col='Year', usecols=lambda x: x not in skipcols)
#this file path navigates from current working directory back up 2 folders, and over to the csv
#convert micron to cm
wafer_thick_cm = wafer_thickness/10000 # microns in a cm
#print(wafer_thick_cm)

# In[17]:

#There are missing data, so we will interpolate linearly for missing years
wafer_thick_cm = wafer_thick_cm.interpolate(method='linear',axis=0)
#print(wafer_thick_cm)
plt.plot(wafer_thick_cm, label='Wafer Thickness (cm)')
plt.title('Wafer Thickness (cm)')

# Now multiply the thickness of the cell by the area of the cell to get a cell volume for each year

# In[18]:

#First, remove 1990 through 1994, to match the size of the cell area df
wafer_thick_cm_1995 = wafer_thick_cm.loc[wafer_thick_cm.index>=1995]
#rename the columns for df.mul operation
wafer_thick_cm_1995.columns = ['avg_cell']
df_cell_volume = df_cellarea_cm2.mul(wafer_thick_cm_1995,'columns')
df_cell_volume.columns = ['cell_volume_cm3']
#print(df_cell_volume)
#plt.plot(df_cell_volume, label='average cell volume (cm^3)')
#plt.legend()

# Now we have the volume of the cell in cm^3 for each year, we can bring in the density of Silicon to get a mass of Silicon per cell for each year.

# In[19]:

df_Simass_percell = df_cell_volume.mul(density_si)
df_Simass_percell.columns = ['Si_gpercell']
print(df_Simass_percell)
df_Simass_percell.to_csv(cwd+'/../../PV_ICE/baselines/SupportingMaterial/output_si_g_per_cell.csv', index=True)
plt.plot(df_Simass_percell, label='Mass Si per cell (g/cell)')
#plt.legend()
plt.title('Mass Silicon per cell annually')
plt.xlabel('Year')
plt.ylabel('Silicon (grams/cell)')

# ## g Si per m^2 of module

# Now take the above mass of silicon per cell and multiply it by the factor of number of cells per m^2 of module

# In[20]:

df_Simass_percell.columns = df_cellperm2.columns = ['Si_g'] #rename to a common name
df_Simass_perm2 = df_Simass_percell.mul(df_cellperm2, 'columns') #multiply
#print(df_Simass_perm2)
#print out to a csv
df_Simass_perm2.to_csv(cwd+'/../../PV_ICE/baselines/SupportingMaterial/output_si_g_per_m2.csv', index=True)

#make a pretty plot
plt.plot(df_Simass_perm2, label='Silicon g/m^2 of module')
plt.legend()
plt.title('Mass of Si per module area')
plt.xlabel('Year')
plt.ylabel('Silicon (grams/m^2)')

# For post-2030, the mass per m^2 of silicon was held constant through 2050 due to the uncertainty about future technology trends. For example, there are at least 3 different cell sizes which are vying for becoming the next mainstream technology, the move to all bifiacial might affect silicon use differently, and the half-cut and smaller cell technologies will also have an affect. Therefore, we have held it constant from 2030 onward, and this assumption can be modified by the user.

# In[21]:

#understanding what influences the changes in Si mass/module m^2
plt.plot(df_Simass_perm2, label='Silicon g/m^2')
plt.plot(wafer_thick_cm, label='Wafer Thickness (cm)')
plt.plot(market_average_cell_dims, label='annual avg cell dims in mm')
plt.plot(scaled_marketshares['monoSi'],label='MonoSi Marketshare')

plt.legend()
plt.title('Mass of Si per module area by constituent parts')
plt.xlim([1995,2030])
plt.xlabel('Year')

max_stars_repo_path: docs/tutorials/baseline development documentation/(baseline development) Silicon per m2.py


text:
import random

import numpy as np
import torch
import torch.jit as jit
import torch.nn as nn
import torch.nn.functional as F

from habitat_baselines.rl.aux_losses import AuxLosses
from habitat_baselines.rl.models.transformer_xl import GELU, Swish, TransformerXL

class TransformerStateEncoder(nn.Module):
    def __init__(
        self,
        input_size: int,
        n_layer: int = 6,
        d_head: int = 64,
        d_model: int = 256,
        d_inner: int = 256,
        layer_drop: float = 0.0,
        dropout: float = 0.0,
        dropatt: float = 0.0,
        self_sup: bool = False,
        max_self_sup_K: int = 30,
    ):

        super().__init__()
        assert input_size == d_model
        n_head = d_model // d_head

        self.transformer = jit.script(
            TransformerXL(
                n_layer=n_layer,
                n_head=n_head,
                d_model=d_model,
                d_head=d_head,
                d_inner=d_inner,
                dropout=dropout,
                dropatt=dropatt,
                layer_drop=layer_drop,
            )
        )

        self.max_self_sup_K = max_self_sup_K
        self.self_sup = self_sup
        self.classifier = nn.Sequential(
            nn.Linear(d_model * 2, d_model // 2, bias=False),
            nn.LayerNorm(d_model // 2),
            nn.ReLU(True),
            nn.Linear(d_model // 2, 1),
        )

        self.layer_init()

    def layer_init(self):
        for name, param in self.transformer.named_parameters():
            if "weight" in name and len(param.size()) >= 2:
                nn.init.orthogonal_(param)
            elif "bias" in name:
                nn.init.constant_(param, 0)

    def single_forward(self, x, context, mems, masks):
        r"""Forward for a non-sequence input
        """
        mems = mems * masks.view(1, -1, 1, 1)
        x, mems, _ = self.transformer(x.unsqueeze(0), context.unsqueeze(0), mems)
        x = x.squeeze(0)

        return x, mems

    def seq_forward(self, x, context, mems, masks):
        r"""Forward for a sequence of length T
        Args:
            x: (T, N, -1) Tensor that has been flattened to (T * N, -1)
            hidden_states: The starting hidden state.
            masks: The masks to be applied to hidden state at every timestep.
                A (T, N) tensor flatten to (T * N)
        """
        # x is a (T, N, -1) tensor flattened to (T * N, -1)
        n = mems.size(1)
        t = int(x.size(0) / n)

        # unflatten
        x = x.view(t, n, x.size(1))
        context = context.view(t, n, context.size(1))
        masks = masks.view(t, n)

        if self.self_sup:
            ep_lens = []
            for i in range(n):
                last_zero = 0
                has_zeros = (
                    (masks[1:-1, i] == 0.0).nonzero().squeeze(-1).cpu().unbind(0)
                )
                for z in has_zeros:
                    z = z.item() + 1
                    ep_lens.append(z - last_zero)
                    last_zero = z

                ep_lens.append(t - last_zero)

            k = random.randint(
                1,
                max(min(self.max_self_sup_K, int(0.8 * np.mean(np.array(ep_lens)))), 2),
            )
        else:
            k = None

        content, mems, query = self.transformer.transformer_seq_forward(
            x, context, mems, masks, two_stream_k=k
        )

        if self.self_sup:
            positives = x
            negatives = []
            for _ in range(3):
                negative_inds = torch.randperm(t * n, device=x.device)
                negatives.append(
                    torch.gather(
                        x.view(t * n, -1),
                        dim=0,
                        index=negative_inds.view(t * n, 1).expand(t * n, x.size(-1)),
                    ).view(t, n, -1)
                )

            negatives = torch.stack(negatives, dim=-1)

            positives = torch.einsum("...i, ...i -> ...", positives, query)
            negatives = torch.einsum("...ik, ...i -> ...k", negatives, query)
            cpc_logits = torch.stack([positives.unsqueeze(-1), negatives], dim=-1)

            valid_modeling_queries = torch.ones(
                t, n, device=query.device, dtype=torch.bool
            )
            valid_modeling_queries[0:k] = 0
            for i in range(n):
                has_zeros_batch = (
                    (masks[:, i] == 0.0).nonzero().squeeze(-1).cpu().unbind(0)
                )
                for z in has_zeros_batch:
                    valid_modeling_queries[n : n + k, i] = 0

            cpc_loss = torch.masked_select(
                F.cross_entropy(
                    cpc_logits,
                    torch.zeros(t, n, dtype=torch.long, device=cpc_logits.device),
                    reduction="none",
                ),
                valid_modeling_queries,
            ).mean()

            AuxLosses.register_loss("CPC|A", cpc_loss, 0.2)

        content = content.view(t * n, -1)

        return content, mems

    def forward(self, x, context, mems, masks):
        if x.size(0) == mems.size(1):
            return self.single_forward(x, context, mems, masks)
        else:
            return self.seq_forward(x, context, mems, masks)

    def initial_hidden(self, bsz):
        return self.transformer.init_mems(bsz)

    @property
    def num_recurrent_layers(self):
        return 1

max_stars_repo_name: NicoGrande/habitat-api

max_issues_repo_name: NicoGrande/habitat-api

max_forks_repo_name: NicoGrande/habitat-api


text:
"""
Utility functions for CSD methods.
"""

import autograd.numpy as np

def normalize(x):
    return x/np.max(np.abs(x), axis=(0, 1))

def sort_grid(x):
    xsrt = x[x[:, 1].argsort()]
    xsrt = xsrt[xsrt[:, 0].argsort(kind='mergesort')]
    return xsrt

def expand_grid(x1,x2):
    """
    Creates (len(x1)*len(x2), 2) array of points from two vectors.
    :param x1: vector 1, (nx1, 1)
    :param x2: vector 2, (nx2, 1)
    :return: (nx1*nx2, 2) points
    """
    lc = [(a, b) for a in x1 for b in x2]
    return np.squeeze(np.array(lc))

def reduce_grid(x):
    """
    Undoes expand_grid to take (nx, 2) array to two vectors containing unique values of each col.
    :param x: (nx, 2) points
    :return: x1, x2 each vectors
    """
    x1 = np.sort(np.unique(x[:,0]))
    x2 = np.sort(np.unique(x[:,1]))
    return x1, x2

def mykron(A, B):
    """
    Efficient Kronecker product.
    """
    a1, a2 = A.shape
    b1, b2 = B.shape
    C = np.reshape(np.expand_dims(A, (1, 3)) * np.expand_dims(B, (0, 2)), (a1*b1, a2*b2))
    return C

def comp_eig_D(Ks, Kt, sig2n):
    """
    Computes eigvecs and diagonal D for inversion of kron(Ks, Kt) + sig2n * I
    :param Ks: spatial covariance
    :param Kt: temporal covariance
    :param sig2n: noise variance
    :return: eigvec(Ks), eigvec(Kt), Dvec
    """
    nx = Ks.shape[0]
    nt = Kt.shape[0]
    if np.isscalar(sig2n):
        sig2n_vec = sig2n*np.ones(nx*nt)
    else:
        sig2n_vec = np.repeat(sig2n, nt) #sig2n can be nx dimension
    evals_t, evec_t = np.linalg.eigh(Kt)
    evals_s, evec_s = np.linalg.eigh(Ks)
    #import scipy.linalg
    #evals_t, evec_t = scipy.linalg.eigh(Kt)
    #evals_s, evec_s = scipy.linalg.eigh(Ks)
    Dvec = np.repeat(evals_s, nt) * np.tile(evals_t, nx) + sig2n_vec
    return evec_s, evec_t, Dvec

max_stars_repo_name: natalieklein/gpcsd

max_issues_repo_name: natalieklein/gpcsd

max_forks_repo_name: natalieklein/gpcsd

reason: import scipy


text:
#!/usr/bin/python()
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import os
import numpy as np
BASE_DIR = os.path.join(os.path.dirname(__file__),'../../../data/Stanford3dDataset_v1.2_Aligned_Version')

object_dict = {
            'clutter':   0,
            'ceiling':   1,
            'floor':     2,
            'wall':      3,
            'beam':      4,
            'column':    5,
            'door':      6,
            'window':    7,
            'table':     8,
            'chair':     9,
            'sofa':     10,
            'bookcase': 11,
            'board':    12}

path_Dir_Areas =  os.listdir(BASE_DIR)

for Area in path_Dir_Areas:
    path_Dir_Rooms = os.listdir(os.path.join(BASE_DIR,Area))
    for Room in path_Dir_Rooms:
        xyz_Room = np.zeros((1,6))
        label_Room = np.zeros((1,1))
        path_Annotations = os.path.join(BASE_DIR,Area,Room,"Annotations")
        print(path_Annotations)
        # make store directories
        path_prepare_label = os.path.join("../../../data/S3DIS/prepare_label_rgb",Area,Room)
        if not os.path.exists(path_prepare_label):
            os.makedirs(path_prepare_label)
        #############################
        path_objects = os.listdir(path_Annotations)
        for Object in path_objects:
            if object_dict.has_key(Object.split("_",1)[0]):
                print(Object.split("_",1)[0] + " value:" ,object_dict[Object.split("_",1)[0]])
                xyz_object = np.loadtxt(os.path.join(path_Annotations,Object))[:,:]#(N,6)
                label_object = np.tile([object_dict[Object.split("_",1)[0]]],(xyz_object.shape[0],1))#(N,1)
            else:
                continue

            xyz_Room = np.vstack((xyz_Room,xyz_object))
            label_Room = np.vstack((label_Room,label_object))

        xyz_Room = np.delete(xyz_Room,[0],0)
        label_Room = np.delete(label_Room,[0],0)

        np.save(path_prepare_label+"/xyzrgb.npy",xyz_Room)
        np.save(path_prepare_label+"/label.npy",label_Room)


text:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import requests
import math

path = './server/'

urlConfirmed = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv'
urlDeaths = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv'
urlCountries = 'https://restcountries.eu/rest/v2/all'
urlCovidData = 'https://covid.ourworldindata.org/data/owid-covid-data.json'
super_file = 'super_file.csv'
#urlTemp = 'https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_3.23/crucy.1506241137.v3.23/countries/tmp/crucy.v3.23.1901.2014.{}.tmp.per'
wldMean = 100

df_cc = pd.read_csv(urlConfirmed)
indexesConfirmed = np.arange(2, len(df_cc.columns) - 1)
df_cc.drop(df_cc.columns[indexesConfirmed], axis=1, inplace=True)
df_cc.drop(df_cc.columns[0], axis=1, inplace=True)
df_cc.columns = ['Name','Cases']
dict_cc = {"Name": [], "Cases": []}

name = df_cc.iloc[0,0]
count = 0

for i in range(0, len(df_cc.index)-1):
    n = df_cc.iloc[i,0]
    if not n == name:
    	if name in dict_cc["Name"]:
    		dict_cc["Cases"][dict_cc["Name"].index(name)] += count
    	else:
	    	dict_cc["Name"].append(name)
	    	dict_cc["Cases"].append(count)
    	name = df_cc.iloc[i, 0]
    	if name == 'US':
    		name = 'United States of America'
    	count = 0
    count += df_cc.iloc[i, 1]

df_cc = pd.DataFrame.from_dict(dict_cc)

df_de = pd.read_csv(urlDeaths)
indexesDeaths = np.arange(2, len(df_de.columns) - 1)
df_de.drop(df_de.columns[indexesDeaths], axis=1, inplace=True)
df_de.drop(df_de.columns[0], axis=1, inplace=True)
df_de.columns = ['Name', 'Deaths']

dict_de = {"Name": [], "Deaths": []}

name = df_de.iloc[0,0]
count = 0

for i in range(0, len(df_de.index)-1):
    n = df_de.iloc[i,0]
    if not n == name:
    	if name in df_de["Name"]:
    		dict_de["Deaths"][dict_de["Name"].index(name)] += count
    	else:
	    	dict_de["Name"].append(name)
	    	dict_de["Deaths"].append(count)
    	name = df_de.iloc[i,0]
    	if name == 'US':
    		name = 'United States of America'
    	count = 0
    count += df_de.iloc[i, 1]

df_de = pd.DataFrame.from_dict(dict_de)
df_de.drop(df_de.columns[0], axis=1, inplace=True)

r = requests.get(urlCountries)
arrCountries = r.json()
names = []
populations = []
areas = []
dds = []
lats = []
longs = []
codes = []

for country in arrCountries:
    name = country["name"]
    population = country["population"]
    area = country["area"]
    latlng = country["latlng"]
    code = country["alpha3Code"]

    if name == None or population == None or area == None or latlng == None or code == None:
        continue

    names.append(name)
    populations.append(population)
    areas.append(area)
    dds.append((population/area))
    lats.append(latlng[0])
    longs.append(latlng[1])
    codes.append(code)

dictCountries = {"Name": names, "Code": codes, "Latitude":lats, "Longitude":longs, "Population": populations, "Area": areas, "Demographic_Density": dds}
df_co = pd.DataFrame.from_dict(dictCountries)

df_super = pd.read_csv(super_file)
dict_super = {"Code": []}
code = df_super.iloc[0, 1]
dict_super["Code"].append(code)

for l in range(0, len(df_super.index)):
    series_name = df_super.iloc[l, 2]
    if not code == df_super.iloc[l, 1]:
        code = df_super.iloc[l, 1]

        if not str(code) == 'nan':
            dict_super["Code"].append(code)

    if not isinstance(series_name, str) or not isinstance(code, str):
        continue

    for c in range(len(df_super.columns)-1, -1, -1):
        value = df_super.iloc[l, c]
        if (isinstance(value, int) and math.isnan(value)):
            continue
        try:
            value = float(value)
        except:
            value = df_super.iloc[l, c]
            if not value == '..':
                if series_name in dict_super:
                    dict_super[series_name].append(np.nan)
                else:
                    dict_super[series_name] = [np.nan]
                break
            continue

        if series_name in dict_super:
            dict_super[series_name].append(value)
        else:
            dict_super[series_name] = [value]

        break

df_super = pd.DataFrame.from_dict(dict_super)
arr = ['Population ages 65 and above (% of total population)', 'Urban population (% of total population)', 'Demographic Density', 'Hospital beds (per 1,000 people)', 'Code']
df_super.drop(df_super.loc[:, ~df_super.columns.isin(arr)], axis=1, inplace=True)
#df_super.drop(['Investment in water and sanitation with private participation (current US$)', 'Community health workers (per 1,000 people)', 'Adequacy of social safety net programs (% of total welfare of beneficiary households)', 'Adequacy of social insurance programs (% of total welfare of beneficiary households)', 'Railways, passengers carried (million passenger-km)', 'Completeness of death registration with cause-of-death information (%)', 'Population in urban agglomerations of more than 1 million (% of total population)', 'Electric power consumption (kWh per capita)', 'Smoking prevalence, total (ages 15+)' ,'Average precipitation in depth (mm per year)'], axis=1, inplace=True)
df_super.drop([len(df_super.index)-1], inplace=True)

r = requests.get(urlCovidData)
dict_json_cpm = r.json()
dict_cpm = {"Code": [], "Deaths_per_mil": []}
for key in dict_json_cpm.keys():
    cpm = dict_json_cpm[key][len(dict_json_cpm[key])-1]
    try:
        cpm = cpm['total_deaths_per_million']
    except:
        cpm = cpm['total_deaths']/(cpm['population']/10**6)
    dict_cpm["Code"].append(key)
    if cpm > wldMean:
        dict_cpm["Deaths_per_mil"].append('Above')
    else:
        dict_cpm["Deaths_per_mil"].append('Below')
df_cpm = pd.DataFrame.from_dict(dict_cpm)

df_final = df_cc
df_final['Deaths'] = df_de['Deaths']
df_final = df_final.join(df_co.set_index('Name'), on='Name')
df_final = df_final.join(df_super.set_index('Code'), on='Code')
df_final = df_final.join(df_cpm.set_index('Code'), on='Code')
df_final.set_index('Name', inplace=True)

df_final = df_final.dropna()

df_final.to_csv(path+'final.csv')


text:
"""
Credits:
Copyright (c) 2017-2019 Matej Aleksandrov, Matej Batič, Andrej Burja, Eva Erzin (Sinergise)
Copyright (c) 2017-2019 Grega Milčinski, Matic Lubej, Devis Peresutti, Jernej Puc, Tomislav Slijepčević (Sinergise)
Copyright (c) 2017-2019 Blaž Sovdat, Jovan Višnjić, Anže Zupanc, Lojze Žust (Sinergise)

This source code is licensed under the MIT license found in the LICENSE
file in the root directory of this source tree.
"""

import pytest
import unittest
import os
import numpy as np

from eolearn.mask import AddCloudMaskTask, get_s2_pixel_cloud_detector
from eolearn.core import EOPatch, FeatureType

class TestAddSentinelHubCloudMaskTask(unittest.TestCase):

    TEST_PATCH_FILENAME = os.path.join(os.path.dirname(os.path.realpath(__file__)), '../../../example_data',
                                       'TestEOPatch')

    FEATURES_TO_LOAD = [
        (FeatureType.DATA, 'BANDS-S2-L1C', 'CLP'),
        (FeatureType.MASK, 'CLM', 'IS_DATA'),
        (FeatureType.LABEL, 'IS_CLOUDLESS'),
        FeatureType.META_INFO,
        FeatureType.TIMESTAMP,
        FeatureType.BBOX
    ]

    @classmethod
    def setUpClass(cls):
        cls.eop = EOPatch.load(cls.TEST_PATCH_FILENAME, features=cls.FEATURES_TO_LOAD)
        cls.eop.rename_feature(FeatureType.DATA, 'BANDS-S2-L1C', 'ALL_DATA')

    def test_raises_errors(self):
        classifier = get_s2_pixel_cloud_detector(all_bands=True)
        add_cm = AddCloudMaskTask(classifier, 'bands', cmask_feature='clm')
        self.assertRaises(ValueError, add_cm, self.eop)

    def _check_shape(self, output, data):
        """Checks that shape of data and output match."""

        self.assertTrue(output.ndim == 4)
        self.assertTrue(output.shape[:-1] == data.shape[:-1])
        self.assertTrue(output.shape[-1] == 1)

    def test_cloud_coverage(self):
        classifier = get_s2_pixel_cloud_detector(all_bands=True)
        # Classifier is run on same resolution as data array
        add_cm = AddCloudMaskTask(classifier, 'ALL_DATA', cmask_feature='CLM_TEST', cprobs_feature='CLP_TEST')
        eop_clm = add_cm(self.eop)

        # Check shape and type
        self._check_shape(eop_clm.mask['CLM_TEST'], eop_clm.data['ALL_DATA'])
        self._check_shape(eop_clm.data['CLP_TEST'], eop_clm.data['ALL_DATA'])
        self.assertTrue(eop_clm.mask['CLM_TEST'].dtype == np.bool)
        self.assertTrue(eop_clm.data['CLP_TEST'].dtype == np.float32)

        # Compare mean cloud coverage with provided reference
        mean_clm_provided = np.mean(eop_clm.mask['CLM'])
        mean_clp_provided = np.mean(eop_clm.data['CLP'])
        self.assertAlmostEqual(np.mean(eop_clm.mask['CLM_TEST']), mean_clm_provided, places=1)
        self.assertAlmostEqual(np.mean(eop_clm.data['CLP_TEST']), mean_clp_provided, places=2)

        # Classifier is run on downscaled version of data array
        add_cm = AddCloudMaskTask(classifier, 'ALL_DATA', cm_size_y='20m', cm_size_x='20m',
                                  cmask_feature='CLM_TEST', cprobs_feature='CLP_TEST')
        eop_clm = add_cm(self.eop)

        # Check shape and type
        self._check_shape(eop_clm.mask['CLM_TEST'], eop_clm.data['ALL_DATA'])
        self._check_shape(eop_clm.data['CLP_TEST'], eop_clm.data['ALL_DATA'])
        self.assertTrue(eop_clm.mask['CLM_TEST'].dtype == np.bool)
        self.assertTrue(eop_clm.data['CLP_TEST'].dtype == np.float32)

        # Compare mean cloud coverage with provided reference
        mean_clm_provided = np.mean(eop_clm.mask['CLM'])
        mean_clp_provided = np.mean(eop_clm.data['CLP'])
        self.assertAlmostEqual(np.mean(eop_clm.mask['CLM_TEST']), mean_clm_provided, places=2)
        self.assertAlmostEqual(np.mean(eop_clm.data['CLP_TEST']), mean_clp_provided, places=2)

        # Check if same times are flagged as cloudless
        cloudless = np.mean(eop_clm.mask['CLM_TEST'], axis=(1,2,3)) == 0
        self.assertTrue(np.all(cloudless == eop_clm.label['IS_CLOUDLESS'][:,0]))

    @pytest.mark.xfail
    def test_wms_request(self):
        classifier = get_s2_pixel_cloud_detector(all_bands=False)
        # Classifier is run on new request of data array
        add_cm = AddCloudMaskTask(classifier, 'BANDS-S2-L1C', cm_size_x='20m', cm_size_y='20m',
                                  cmask_feature='CLM_TEST', cprobs_feature='CLP_TEST')
        eop_clm = add_cm(self.eop)

        # Check shape and type
        self._check_shape(eop_clm.mask['CLM_TEST'], eop_clm.data['ALL_DATA'])
        self._check_shape(eop_clm.data['CLP_TEST'], eop_clm.data['ALL_DATA'])
        self.assertTrue(eop_clm.mask['CLM_TEST'].dtype == np.bool)
        self.assertTrue(eop_clm.data['CLP_TEST'].dtype == np.float32)

        # Compare mean cloud coverage with provided reference
        mean_clm_provided = np.mean(eop_clm.mask['CLM'])
        mean_clp_provided = np.mean(eop_clm.data['CLP'])
        self.assertAlmostEqual(np.mean(eop_clm.mask['CLM_TEST']), mean_clm_provided, places=2)
        self.assertAlmostEqual(np.mean(eop_clm.data['CLP_TEST']), mean_clp_provided, places=2)

if __name__ == '__main__':
    unittest.main()


text:
import numpy as np
import theano
import logging

from numpy.testing import assert_allclose
from theano import tensor, function

from blocks import initialization
from blocks.bricks import Identity, Linear
from blocks.bricks.recurrent import SimpleRecurrent, LSTM
from blocks.initialization import Constant

from blocks_contrib.bricks.recurrent import DelayLine
from blocks_contrib.bricks.recurrent import Unfolder, UnfolderLSTM

logger = logging.getLogger(__name__)

def test_constant_input_lstm():
    x = tensor.matrix('x')
    proto = LSTM(activation=Identity(), dim=1,
                weights_init=Constant(1/4.),
                biases_init=Constant(0.))
    proto.initialize()

    flagger = Linear(input_dim=1, output_dim=1,
                     weights_init=Constant(1./2.),
                     biases_init=Constant(0.))
    flagger.initialize()

    inp2hid = Linear(input_dim=1, output_dim=4,
                     weights_init=Constant(1/4.),
                     biases_init=Constant(0))
    inp2hid.initialize()

    rnn = UnfolderLSTM(proto, flagger)
    rnn.initialize()

    h = inp2hid.apply(x)
    y = rnn.apply(inputs=h, n_steps=10, batch_size=5)

    F = function([x],y)
    X = np.ones((5,1)).astype(theano.config.floatX)

    H = function([x], flagger.apply(x))
    T = H(X)
    #print T

    Y = F(X)
    print Y
    print Y[0].shape

    assert Y[0].shape == (4,5,1)
    #target = np.cumsum(np.ones((6,1,1)),axis=0)
    #assert_allclose(Y, target)

def test_constant_input():
    x = tensor.matrix('x')
    proto = SimpleRecurrent(activation=Identity(), dim=1,
            weights_init=initialization.Identity(1.))
    proto.initialize()

    flagger = Linear(input_dim=1, output_dim=1,
                     weights_init=Constant(1/10.),
                     biases_init=Constant(0.))
    flagger.initialize()

    rnn = Unfolder(proto, flagger)
    rnn.initialize()

    y = rnn.apply(inputs=x, n_steps=10, batch_size=1)

    F = function([x],y)
    X = np.ones((1,1)).astype(theano.config.floatX)

    H = function([x], flagger.apply(x))
    T = H(X)
    print T

    Y = F(X)
    print Y

    target = np.cumsum(np.ones((5,1,1)),axis=0)
    assert_allclose(Y[0], target)

def test_delay_line():
    x = tensor.tensor3('x')
    input_dim = 1
    batch_size = 1
    memory_size = 3
    time_len = 4
    delay_line = DelayLine(input_dim, memory_size,
                           weights_init=Constant(1.))
    delay_line.initialize()
    y = delay_line.apply(x, iterate=True)
    func = function([x], y)

    x_val = np.zeros((time_len, batch_size, input_dim))
    val = np.arange(4)
    #val = np.tile(x_val[np.newaxis], batch_size).T
    x_val[:,0,0] = val.astype(theano.config.floatX)
    x_val = x_val.astype(theano.config.floatX)
    y_val = func(x_val).astype(theano.config.floatX)
    print y_val
    assert y_val.shape == (4,1,3)

if __name__ == '__main__':
    #test_delay_line()
    test_constant_input()
    test_constant_input_lstm()

reason: import numpy,from numpy,import theano,from theano


text:
#
# EOBTidalExternal
#
#  Provides interface to Bernuzzi EOB model
#
# REQUIRES
#    - matlab
#    - run_goeob_standalone.sh  # created by EOB code interface
#
# PROVIDES
#    - hlmoft
#
# COMPARE TO
#    NRWaveformCatalogManager3   : very similar interface

debug_output = False

import numpy as np
import os
import sys
import shutil
import lalsimulation as lalsim
import lal

from scipy.interpolate import interp1d

import pickle

from .. import lalsimutils

rosUseArchivedWaveforms = True

rosDebug = False
#dirBaseFiles =os.environ["HOME"] + "/unixhome/Projects/LIGO-ILE-Applications/ILE-Tides/MatlabCodePolished"
dirBaseFiles =os.environ["EOB_BASE"]
dirBaseFilesArchive =os.environ["EOB_ARCHIVE"]
dirBaseMatlab=os.environ["MATLAB_BASE"]

#default_interpolation_kind = 'quadratic'  # spline interpolation   # very slow!
default_interpolation_kind = 'linear'  # spline interpolation   # very slow!

internal_ModesAvailable = [(2,2), (2,1), (2,-2), (2,-1), (3,3), (3,2), (3,1), (3,-3), (3,-2), (3,-1)]

MsunInSec = lal.MSUN_SI*lal.G_SI/lal.C_SI**3

def myzero(arg):
    return 0
def RangeWrap1d(bound, val,fn):
    return fn   # IDEALLY not necessary with modern interp1d

def ModeToString(pair):
    return str(pair[0])+str(pair[1])   # this is only used for POSITIVE l,m (single digit)

class WaveformModeCatalog:
    """
    Class containing EOB tidal harmonics,  both in dimensionless and dimensional form
    """

    def __init__(self, P,  lmax=2,
                 align_at_peak_l2_m2_emission=True, mode_list_to_load=[],build_fourier_time_window=1000,clean_with_taper=True,use_internal_interpolation_deltaT=None,build_strain_and_conserve_memory=False,reference_phase_at_peak=None,fix_phase_zero_at_coordinate=False):
        self.P  = P
        self.quantity = "h"
        self.fOrbitLower =0.    #  Used to clean results.  Based on the phase of the 22 mode
        self.fMinMode ={}
        # Mode storage convention
        #   - event time in first element
        #     - t=0 MAY be at peak, if opts.align_at_peak_l2m2_emission.  Cannot promise...but raw samples show the valid range
        #     - t in seconds
        #   - h*r/M in second element
        self.waveform_modes = {}
        self.waveform_modes_uniform_in_time={}
        self.waveform_modes_nonuniform_smallest_timestep = {}
        self.waveform_modes_nonuniform_largest_timestep = {}
        self.waveform_modes[(2,2)] = []
        self.waveform_modes_complex = {}
        self.waveform_modes_complex_interpolated = {}
        self.waveform_modes_complex_interpolated_amplitude = {}
        self.waveform_modes_complex_interpolated_phase = {}
        self.waveform_modes_complex_padded = {}
        self.waveform_modes_fourier = {}

        # Require spins zero for now
        if any([P.s1x,P.s1y,P.s1z,P.s2x,P.s2y,P.s2z]):
            print(" FAILURE: Tidal code assumes a nonspinning approximant for now")

        # Run the external script with the necessary parameters
        # Assume environment set correctly
        # Delete files from previous case"
        m1InMsun = P.m1/lal.MSUN_SI
        m2InMsun = P.m2/lal.MSUN_SI
        m1InMsun, m2InMsun = reversed(sorted([m1InMsun, m2InMsun]))   # FORCE m1 > m2

        # Convert lambda's to kappa
        # Note we are passed lambda1, lambda2.  We convert to (l=2) kappaA, kappaB
        #  - Generate lambdatilde
        #  - generate lambda2
        LambdaTilde, DeltaLambdaTilde = lalsimutils.tidal_lambda_tilde(m1InMsun, m2InMsun, P.lambda1,P.lambda2)
        kappaA2, kappaB2 =  lalsimutils.barlamdel_to_kappal( m1InMsun/m2InMsun, LambdaTilde,2)

        fname_base = "stored_"+"_".join([str(m1InMsun),str(m2InMsun),str(kappaA2),str(kappaB2),str(P.fmin)]) +".dir";
        print("  Saving to file (beware collisions!) ", fname_base)

        retrieve_directory = ''
        cwd = os.getcwd();
        if rosUseArchivedWaveforms and os.path.exists(dirBaseFilesArchive+"/"+fname_base) and os.path.exists(dirBaseFilesArchive+"/"+fname_base+"/test_h22.dat"):
            retrieve_directory = dirBaseFilesArchive+"/"+fname_base
            print(" Attempting to use archived waveform data  in ", retrieve_directory)
        else:
            retrieve_directory = dirBaseFilesArchive+"/"+fname_base + "/"
            cmd = dirBaseFiles+"/run_goeob_standalone.sh " + dirBaseMatlab + " " +str(m1InMsun) + " "+str(m2InMsun) + " " + str(kappaA2) + " " + str(kappaB2) + " " + str(P.fmin) + " " + retrieve_directory
            # Create directory
            if not os.path.exists(retrieve_directory):
                print(" Making directory to archive this run ... ", retrieve_directory)
                os.makedirs(retrieve_directory)
            print(" Generating tidal EOB with ", cmd)
            os.chdir(dirBaseFiles); os.system(cmd);
            # retrieve_directory = dirBaseFiles
            # if rosUseArchivedWaveforms:
            #     retrieve_directory = dirBaseFilesArchive+"/"+fname_base
            #     # Create directory
            #     if not os.path.exists(retrieve_directory):
            #         print " Making directory to archive this run ... ", retrieve_directory
            #         os.makedirs(retrieve_directory)
            #     # Copy
            #     print " Copying files into storage directory ", retrieve_directory
            #     shutil.copyfile(dirBaseFiles+"/test_h22.dat", retrieve_directory+"/test_h22.dat")
            #     shutil.copyfile(dirBaseFiles+"/test_h21.dat", retrieve_directory+"/test_h21.dat")
            #     shutil.copyfile(dirBaseFiles+"/test_h33.dat", retrieve_directory+"/test_h33.dat")
            #     shutil.copyfile(dirBaseFiles+"/test_h32.dat", retrieve_directory+"/test_h32.dat")
            #     shutil.copyfile(dirBaseFiles+"/test_h31.dat", retrieve_directory+"/test_h31.dat")

        # First loop: Create all the basic mode data
        # This should ALREADY BE IN PHYSICAL TIME UNITS but have UNPHYSICAL distance scales
        nu = lalsimutils.symRatio(P.m1,P.m2)
        delta = (m1InMsun- m2InMsun)/(m1InMsun+m2InMsun)
        for mode in internal_ModesAvailable:
            if mode[0]<= lmax:
                # Note all modes are REFLECTION-SYMMETRIC
                self.waveform_modes[mode] = np.loadtxt(retrieve_directory+"/test_h"+ModeToString( (mode[0], np.abs(mode[1])))+".dat")
                self.waveform_modes_nonuniform_smallest_timestep[mode] = self.waveform_modes[mode][1,0]-self.waveform_modes[mode][0,0]  # uniform in time
                self.waveform_modes_nonuniform_largest_timestep[mode] = self.waveform_modes[mode][1,0]-self.waveform_modes[mode][0,0]  # uniform in time
                self.waveform_modes_uniform_in_time[mode] =True

                # Add factor of 'nu' that was missing
                if mode[1] %2 ==0 :
                    self.waveform_modes[mode][:,1]*= nu
                    self.waveform_modes[mode][:,2]*= nu
                else:
                    self.waveform_modes[mode][:,1]*= nu*delta
                    self.waveform_modes[mode][:,2]*= nu*delta

                # Create, if symmetric
                if mode[1]<0:  # (-1)^l conjugate
                      self.waveform_modes[mode][:,1] *= (-1)**mode[0]
                      self.waveform_modes[mode][:,2] *= -1.*(-1)**mode[0]

                if clean_with_taper:
                    # chop off the first 1 seconds.
                    # use *physical time* for tapering, because it is nonuniform
                    ncut = np.argmax( - (1-self.waveform_modes[mode][:,0])**2)
                    tcut = float(self.waveform_modes[mode][ncut,0])
                    vectaper = 0.5-0.5*np.cos(np.pi* self.waveform_modes[mode][:ncut,0]/tcut)
                    self.waveform_modes[mode][:ncut,1]*= vectaper
                    self.waveform_modes[mode][:ncut,2]*= vectaper

                self.waveform_modes_complex[mode] = np.array([self.waveform_modes[mode][:,0], self.waveform_modes[mode][:,1]+1.j*self.waveform_modes[mode][:,2]]).T

                if rosDebug:
                    print("  Loaded ",  mode, "; length = ", len(self.waveform_modes[mode]), " processed interval = ", self.waveform_modes[mode][0,0], self.waveform_modes[mode][-1,0], " sampling interval at start ",  self.waveform_modes[mode][1,0]-self.waveform_modes[mode][0,0], " and smallest interval = ", self.waveform_modes_nonuniform_smallest_timestep[mode])

        nmax_orig = np.argmax(np.abs(self.waveform_modes_complex[(2,2)][:,1]))
        tmax_orig = float(self.waveform_modes_complex[(2,2)][nmax_orig,0])

        for mode in internal_ModesAvailable:
           if mode[0]<=lmax:
                # shift all times, if necessary
                if align_at_peak_l2_m2_emission:
                    self.waveform_modes[mode][:,0] += -tmax_orig
                    self.waveform_modes_complex[mode][:,0] += -tmax_orig
                # Either store in a different field OR reformat
                # Create an interpolating function for the complex amplitude and phase
                datAmp = np.abs(self.waveform_modes_complex[mode][:,1])
                datPhase = lalsimutils.unwind_phase(np.angle(self.waveform_modes_complex[mode][:,1]))    #needs to be unwound to be continuous
                datT = self.waveform_modes_complex[mode][:,0]
                # Compute starting frequency for this mode (used for FFT)
                # Index offset corresponding to 10 M
                dt = self.waveform_modes_complex[mode][1,0]-self.waveform_modes_complex[mode][0,0]
                nOffsetForPhase = 5  # ad-hoc offset based on uniform sampling
                nStride = 5
                self.fMinMode[mode] = np.abs((datPhase[nOffsetForPhase+nStride]-datPhase[nOffsetForPhase])/(2*np.pi*(datT[nOffsetForPhase+nStride]-datT[nOffsetForPhase]))) # for FFT. Modified to work with nonuniform sampling
                # WARNING: this is really only reliable for non-junky resolved modes.  We will clean this lower frequency later.
                if mode ==(2,2):
                    self.fOrbitLower  = 0.5*self.fMinMode[mode]

                # Interpolate. (Works even for nonuniform time sampling)
                if True: #not build_strain_and_conserve_memory:
                    # this is gravy - not critical. BUT it is critical for nonuniform
                    self.waveform_modes_complex_interpolated_amplitude[mode] = RangeWrap1d([np.min(self.waveform_modes_complex[mode][:,0]), np.max(self.waveform_modes_complex[mode][:,0])], 0, interp1d( self.waveform_modes_complex[mode][:,0].astype(float), datAmp,kind=default_interpolation_kind, fill_value=0.,bounds_error=False))
                    self.waveform_modes_complex_interpolated_phase[mode] = RangeWrap1d([np.min(self.waveform_modes_complex[mode][:,0]), np.max(self.waveform_modes_complex[mode][:,0])], 0,  interp1d(self.waveform_modes_complex[mode][:,0].astype(float), datPhase,kind=default_interpolation_kind, fill_value=0.,bounds_error=False))

        print(" Restoring current working directory... ",cwd)
        os.chdir(cwd);

    def complex_hoft(self,  force_T=False, deltaT=1./16384, time_over_M_zero=0.,sgn=-1):
        hlmT = self.hlmoft( force_T, deltaT,time_over_M_zero)
        npts = hlmT[(2,2)].data.length
        wfmTS = lal.CreateCOMPLEX16TimeSeries("Psi4", lal.LIGOTimeGPS(0.), 0., deltaT, lalsimutils.lsu_DimensionlessUnit, npts)
        wfmTS.data.data[:] = 0   # SHOULD NOT BE NECESARY, but the creation operator doesn't robustly clean memory
        wfmTS.epoch = hlmT[(2,2)].epoch
        for mode in hlmT.keys():
            # PROBLEM: Be careful with interpretation. The incl and phiref terms are NOT tied to L.
            if rosDebug:
                print(mode, np.max(hlmT[mode].data.data), " running max ",  np.max(np.abs(wfmTS.data.data)))
            wfmTS.data.data += np.exp(-2*sgn*1j*self.P.psi)* hlmT[mode].data.data*lal.SpinWeightedSphericalHarmonic(self.P.incl,-self.P.phiref,-2, int(mode[0]),int(mode[1]))
        return wfmTS
    def complex_hoff(self, force_T=False):
        htC  = self.complex_hoft( force_T=force_T,deltaT= self.P.deltaT)
        TDlen = int(1./self.P.deltaF * 1./self.P.deltaT)
        assert TDlen == htC.data.length
        hf = lal.CreateCOMPLEX16FrequencySeries("Template h(f)",
                                                htC.epoch, htC.f0, 1./htC.deltaT/htC.data.length, lalsimutils.lsu_HertzUnit,
                                                htC.data.length)
        fwdplan=lal.CreateForwardCOMPLEX16FFTPlan(htC.data.length,0)
        lal.COMPLEX16TimeFreqFFT(hf, htC, fwdplan)
        return hf
    def real_hoft(self,Fp=None, Fc=None):
        """
        Returns the real-valued h(t) that would be produced in a single instrument.
        Translates epoch as needed.
        Based on 'hoft' in lalsimutils.py
        """
        # Create complex timessereis
        htC = self.complex_hoft(force_T=1./self.P.deltaF, deltaT= self.P.deltaT)  # note P.tref is NOT used in the low-level code
        TDlen  = htC.data.length
        if rosDebug:
            print("Size sanity check ", TDlen, 1/(self.P.deltaF*self.P.deltaT))
            print(" Raw complex magnitude , ", np.max(htC.data.data))

        # Create working buffers to extract data from it -- wasteful.
        hp = lal.CreateREAL8TimeSeries("h(t)", htC.epoch, 0.,
            self.P.deltaT, lalsimutils.lsu_DimensionlessUnit, TDlen)
        hc = lal.CreateREAL8TimeSeries("h(t)", htC.epoch, 0.,
            self.P.deltaT, lalsimutils.lsu_DimensionlessUnit, TDlen)
        hT = lal.CreateREAL8TimeSeries("h(t)", htC.epoch, 0.,
            self.P.deltaT, lalsimutils.lsu_DimensionlessUnit, TDlen)
        # Copy data components over
        hp.data.data = np.real(htC.data.data)
        hc.data.data = np.imag(htC.data.data)
        # transform as in lalsimutils.hoft
        if Fp!=None and Fc!=None:
            hp.data.data *= Fp
            hc.data.data *= Fc
            hp = lal.AddREAL8TimeSeries(hp, hc)
            hoft = hp
        elif self.P.radec==False:
            fp = lalsimutils.Fplus(self.P.theta, self.P.phi, self.P.psi)
            fc = lalsimutils.Fcross(self.P.theta, self.P.phi, self.P.psi)
            hp.data.data *= fp
            hc.data.data *= fc
            hp.data.data  = lal.AddREAL8TimeSeries(hp, hc)
            hoft = hp
        else:
            # Note epoch must be applied FIRST, to make sure the correct event time is being used to construct the modulation functions
            hp.epoch = hp.epoch + self.P.tref
            hc.epoch = hc.epoch + self.P.tref
            if rosDebug:
                print(" Real h(t) before detector weighting, ", np.max(hp.data.data), np.max(hc.data.data))
            hoft = lalsim.SimDetectorStrainREAL8TimeSeries(hp, hc,    # beware, this MAY alter the series length??
                self.P.phi, self.P.theta, self.P.psi,
                lalsim.DetectorPrefixToLALDetector(str(self.P.detector)))
            hoft = lal.CutREAL8TimeSeries(hoft, 0, hp.data.length)       # force same length as before??
            if rosDebug:
                print("Size before and after detector weighting " , hp.data.length, hoft.data.length)
        if rosDebug:
            print(" Real h_{IFO}(t) generated, pre-taper : max strain =", np.max(hoft.data.data))
        if self.P.taper != lalsimutils.lsu_TAPER_NONE: # Taper if requested
            lalsim.SimInspiralREAL8WaveTaper(hoft.data, self.P.taper)
        if self.P.deltaF is not None:
            TDlen = int(1./self.P.deltaF * 1./self.P.deltaT)
            print("Size sanity check 2 ", int(1./self.P.deltaF * 1./self.P.deltaT), hoft.data.length)
            assert TDlen >= hoft.data.length
            npts = hoft.data.length
            hoft = lal.ResizeREAL8TimeSeries(hoft, 0, TDlen)
            # Zero out the last few data elements -- NOT always reliable for all architectures; SHOULD NOT BE NECESSARY
            hoft.data.data[npts:TDlen] = 0

        if rosDebug:
            print(" Real h_{IFO}(t) generated : max strain =", np.max(hoft.data.data))
        return hoft

    def non_herm_hoff(self):
        """
        Returns the 2-sided h(f) associated with the real-valued h(t) seen in a real instrument.
        Translates epoch as needed.
        Based on 'non_herm_hoff' in lalsimutils.py
        """
        htR = self.real_hoft() # Generate real-valued TD waveform, including detector response
        if self.P.deltaF == None: # h(t) was not zero-padded, so do it now
            TDlen = nextPow2(htR.data.length)
            htR = lal.ResizeREAL8TimeSeries(htR, 0, TDlen)
        else: # Check zero-padding was done to expected length
            TDlen = int(1./self.P.deltaF * 1./self.P.deltaT)
            assert TDlen == htR.data.length
        fwdplan=lal.CreateForwardCOMPLEX16FFTPlan(htR.data.length,0)
        htC = lal.CreateCOMPLEX16TimeSeries("hoft", htR.epoch, htR.f0,
            htR.deltaT, htR.sampleUnits, htR.data.length)
        # copy h(t) into a COMPLEX16 array which happens to be purely real
        htC.data.data[:htR.data.length] = htR.data.data
#        for i in range(htR.data.length):
#            htC.data.data[i] = htR.data.data[i]
        hf = lal.CreateCOMPLEX16FrequencySeries("Template h(f)",
            htR.epoch, htR.f0, 1./htR.deltaT/htR.data.length, lalsimutils.lsu_HertzUnit,
            htR.data.length)
        lal.COMPLEX16TimeFreqFFT(hf, htC, fwdplan)
        return hf

    def estimateFminHz(self):
        return 2*self.fOrbitLower/(MsunInSec*(self.P.m1+self.P.m2)/lal.MSUN_SI)

    def estimateDurationSec(self):
        """
        estimateDuration uses the ACTUAL UNITS IN THE WAVEFORM, which are already in sec
        """
        return np.real(self.waveform_modes_complex[(2,2)][-1,0]-self.waveform_modes_complex[(2,2)][0,0]) # self.deltaToverM*(self.len(self.waveform_modes[(2,2)])

    def hlmoft(self,  force_T=False, deltaT=1./16384, time_over_M_zero=0.,taper_start_time=True):
        """
        hlmoft uses stored interpolated values for hlm(t) generated via the standard cleaning process, scaling them
        to physical units for use in injection code.

        If the time window is sufficiently short, the result is NOT tapered (!!) -- no additional tapering is applied

        The code will ALWAYS have zero padding on the end -- half of the buffer is zero padding!
        This can cause loss of frequency content if you are not careful
        """
        hlmT ={}
        # Define units
        m_total_s = MsunInSec*(self.P.m1+self.P.m2)/lal.MSUN_SI
        distance_s = self.P.dist/lal.C_SI  # insures valid units

        # Create a suitable set of time samples.  Zero pad to 2^n samples.
        # Note waveform is stored in s already
        T_estimated = np.real(self.waveform_modes_complex[(2,2)][-1,0] - self.waveform_modes_complex[(2,2)][0,0])
        npts=0
        n_crit = 0
        if not force_T:
            npts_estimated = int(T_estimated/deltaT)
#            print " Estimated length: ",npts_estimated, T_estimated
            npts = lalsimutils.nextPow2(npts_estimated)
        else:
            npts = int(force_T/deltaT)
            print(" Forcing length T=", force_T, " length ", npts)
        # WARNING: Time range may not cover the necessary time elements.
        # Plan on having a few seconds buffer at the end
        T_buffer_required = npts*deltaT
        print(" EOB internal: Estimated time window (sec) ", T_estimated, " versus buffer duration ", T_buffer_required)
        print(" EOB internal: Requested size vs buffer size",   npts, len(self.waveform_modes_complex[(2,2)]))
        # If the waveform is longer than the buffer, we need to avoid wraparound

        # If the buffer requested is SHORTER than the 2*waveform, work backwards
        # If the buffer requested is LONGER than the waveform, work forwards from the start of all data
        if T_buffer_required/2 > T_estimated:
            tvals =  np.arange(npts)*deltaT + float(self.waveform_modes_complex[(2,2)][0,0])   # start at time t=0 and go forwards (zeros automatically padded by interpolation code)
            t_crit = float( -self.waveform_modes_complex[(2,2)][0,0])
            n_crit = int( t_crit/deltaT) # estiamted peak sample location in the t array, working forward

        else:
            print("  EOB internal: Warning LOSSY conversion to insure half of data is zeros ")
            # Create time samples by walking backwards from the last sample of the waveform, a suitable duration
            # ASSUME we are running in a configuration with align_at_peak_l2m2_emission
            # FIXME: Change this
            tvals = T_buffer_required/2  + (-npts + 1+ np.arange(npts))*deltaT + np.real(self.waveform_modes_complex[(2,2)][-1,0])  # last insures we get some ringdown
            t_crit = T_buffer_required/2 - (np.real(self.waveform_modes_complex[(2,2)][-1,0]))
            n_crit = int(t_crit/deltaT)

        if rosDebug:
            print(" time range being sampled ", [min(tvals),max(tvals)], " corresponding to dimensionless range", [min(tvals)/m_total_s,max(tvals)/m_total_s])
            print(" estimated peak sample at ", n_crit)

        # Loop over all modes in the system
        for mode in self.waveform_modes.keys():
            amp_vals = m_total_s/distance_s * self.waveform_modes_complex_interpolated_amplitude[mode](tvals)  # vectorized interpolation with piecewise
            phase_vals = self.waveform_modes_complex_interpolated_phase[mode]( tvals)
            phase_vals = lalsimutils.unwind_phase(phase_vals)

            if rosDebug:
                print("  Mode ", mode, " physical strain max, indx,", np.max(amp_vals), np.argmax(amp_vals))

            # Copy into a new LIGO time series object
            wfmTS = lal.CreateCOMPLEX16TimeSeries("h", lal.LIGOTimeGPS(0.), 0., deltaT, lalsimutils.lsu_DimensionlessUnit, npts)
            wfmTS.data.data =amp_vals*np.exp(1j*phase_vals)

            # Set the epoch for the time series correctly: should have peak near center of series by construction
            # note all have the same length
            # wfmTS.epoch = -deltaT*wfmTS.data.length/2  # did not work
            #n_crit = np.argmax(wfmTS.data.data)
            #print n_crit*wfmTS.deltaT, wfmTS.epoch   # this should be nearly zero

            # taper the start (1s. Only needed if I do not grab the whole range, because I taper the raw data)
            if taper_start_time:
                tTaper = 1
                nTaper = int(tTaper/deltaT)
                hoft_window = lal.CreateTukeyREAL8Window(nTaper*2, 0.8)
                factorTaper = hoft_window.data.data[0:nTaper]
                wfmTS.data.data[:nTaper]*=factorTaper

            # Store the resulting mode
            hlmT[mode] = wfmTS

        # Set time at peak of 22 mode. This is a hack, but good enough for us
#        n_crit = np.argmax(hlmT[(2,2)].data.data)
        epoch_crit = float(-t_crit) #-deltaT*n_crit
        print(" EOB internal: zero epoch sample location", n_crit, np.argmax(np.abs(hlmT[(2,2)].data.data)))
        for mode in hlmT:
            hlmT[mode].epoch = epoch_crit

        return hlmT

    def hlmoff(self, force_T=False, deltaT=1./16384, time_over_M_zero=0.):
        """
        hlmoff takes fourier transforms of LAL timeseries generated from hlmoft.
        All modes have physical units, appropriate to a physical signal.
        """
        hlmF ={}
        hlmT = self.hlmoft(force_T=force_T,deltaT=deltaT,time_over_M_zero=time_over_M_zero)
        for mode in self.waveform_modes.keys():
            wfmTS=hlmT[mode]
            # Take the fourier transform
            wfmFD = lalsimutils.DataFourier(wfmTS)  # this creates a new f series for *each* call.
            # Store the resulting mode
            hlmF[mode] = wfmFD
        return hlmF

class WaveformMode:
    """
    Class representing a dimensionless timeseries X(\tau=t/M)
    """

###
### FFT syntatic sugar
###
def DataFourierNumpy(wfComplex):    # assume (n,2) size array of [tvals, g(t)]; return fvals, tilde(g).
    # FFT
    # Unroll -- it will save me time later if we are continuous
    T = wfComplex[-1,0] - wfComplex[0,0]
    n = len(wfComplex)
    dt = T/n
    gtilde = np.fft.fft(wfComplex[:,1])*dt
    # Option 1: Raw pairing, not caring about continuity in the frequency array.
#    fvals = 1./T* np.array([ k if  k<=n/2 else k-n for k in np.arange(n)])
    # Option 2: Unroll
    gtilde = np.roll(gtilde, n/2)
    fvals = 1./T*(np.arange(n) - n/2+1)     # consistency with DataFourier: \int dt e(i2\pi f t), note REVERSING
    wfComplexF = np.array([fvals,gtilde[::-1]]).T
    return wfComplexF

def DataInverseFourierNumpy(wfComplex):    # assume (n,2) size array of [fvals, gtilde(f)]; return tvals, g
#    print "NAN check ", wfComplex[np.isnan(wfComplex[:,1]),1]
    df = wfComplex[1,0] - wfComplex[0,0]
    n = len(wfComplex)
    T = 1./df
    dt = T/n
    # Undo the roll, then perform the FFT
    datReversed = wfComplex[:,1][::-1]
    g = np.fft.ifft(np.roll(datReversed, -n/2))*n*df # undo the reverse and roll
    # Assign correct time values.
    #    - Note the zero of time is now centered in the array -- we don't carry a time reference with us.
    tvals = dt*(np.arange(n) -n/2+1)
    wfComplexT = np.array([tvals,g]).T
    return wfComplexT

###
### Mode syntactic sugar
###
def RawGetModePeakTime(wfMode):   # assumed applied to complex data sequence
    nmax = np.argmax(np.abs(wfMode[:,1]))
    return np.real(wfMode[nmax][0])

###
### Plotting syntactic sugar
###

# http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.plot
mode_line_style= {}
mode_line_style[(2,2)] = 'r'
mode_line_style[(2,1)] = 'b'
mode_line_style[(2,-1)] = 'b'
mode_line_style[(2,-2)] = 'r'
mode_line_style[(3,3)] = 'r-'
mode_line_style[(3,-3)] = 'r-'
mode_line_style[(3,1)] = 'b-'
mode_line_style[(3,-1)] = 'b-'

mode_line_style[(2,0)] = 'g'

for l in np.arange(2,5):
    for m in np.arange(-l,l+1):
        if not ((l,m) in mode_line_style.keys()):
            mode_line_style[(l,m)] = 'k'

max_stars_repo_name: spfanning/research-projects-RIT

max_issues_repo_name: spfanning/research-projects-RIT

max_forks_repo_name: spfanning/research-projects-RIT

reason: import numpy,from scipy


text:
import unittest
import nose.tools
import numpy as np

from ..test_utils import generate_linear_structure

from mia.features._orientated_bins import orientated_bins
from mia.features._nonmaximum_suppression import nonmaximum_suppression

class NonMaximumSuppressionTest(unittest.TestCase):

    def test_suppression(self):
        linear_structure = generate_linear_structure(20)
        nbins = 8
        line_strength, line_orinetation = orientated_bins(linear_structure, nbins)
        suppressed = nonmaximum_suppression(line_strength, line_orinetation, nbins)

        nose.tools.assert_equal(np.count_nonzero(suppressed), 60)

    def test_suppression_fails_incorrect_bins(self):
        linear_structure = generate_linear_structure(20)
        nbins = 6
        line_strength, line_orinetation = orientated_bins(linear_structure, 7,
                                                          nbins=nbins)

        args = (line_strength, line_orinetation, nbins)
        nose.tools.assert_raises(ValueError, nonmaximum_suppression, *args)

    def test_different_kernel_size(self):
        linear_structure = generate_linear_structure(20)
        nbins = 8
        line_strength, line_orinetation = orientated_bins(linear_structure, nbins)
        suppressed = nonmaximum_suppression(line_strength, line_orinetation,
                                            nbins, kernel_size=6)

        nose.tools.assert_equal(np.count_nonzero(suppressed), 60)


text:
'''
Date: 2021-08-02 22:38:28
LastEditors: xgy
LastEditTime: 2021-08-15 16:19:23
FilePath: \code\ctpn\src\create_dataset.py
'''
# create train dataset of mindspore
from __future__ import division
import os
import numpy as np
from PIL import Image
from mindspore.mindrecord import FileWriter
from src.config import config

def create_label(train_img_dir, train_txt_dir, prefix=''):
    image_files = []
    image_anno_dict = {}
    img_basenames = []
    for file_name in os.listdir(train_img_dir):
        if 'gif' not in file_name:
            img_basenames.append(os.path.basename(file_name))
    img_names = []
    for item in img_basenames:
        temp1, _ = os.path.splitext(item)
        img_names.append((temp1, item))
    for img, img_basename in img_names:
        image_path = train_img_dir + '/' + img_basename
        annos = []
        file_name = prefix + img + ".txt"
        file_path = os.path.join(train_txt_dir, file_name)
        gt = open(file_path, 'r', encoding='gbk').read().splitlines()
        if not gt:
            continue
        for img_each_label in gt:
            spt = img_each_label.replace(',', '').split(' ')
            if ' ' not in img_each_label:
                spt = img_each_label.split(',')
            annos.append([spt[0], spt[1], spt[2], spt[3]] + [1])
        if annos:
            image_anno_dict[image_path] = np.array(annos)
            image_files.append(image_path)
    return image_files, image_anno_dict

def create_train_dataset(dataset_type):
    image_files = []
    image_anno_dict = {}
    if dataset_type == "pretraining":
        pretrain_image_files, pretrain_anno_dict = create_label(config.hwdb_pretrain_path[0],config.hwdb_pretrain_path[1])
        data_to_mindrecord_byte_image(pretrain_image_files, pretrain_anno_dict,config.pretrain_dataset_path,prefix="ctpn_pretrain.mindrecord", file_num=8)
    elif dataset_type == "finetune":
        finetune_image_files, finetune_anno_dict = create_label(config.hwdb_finetune_path[0],config.hwdb_finetune_path[1])
        data_to_mindrecord_byte_image(finetune_image_files, finetune_anno_dict,config.finetune_dataset_path,prefix="ctpn_finetune.mindrecord", file_num=8)
    elif dataset_type == "test":
        test_image_files, test_anno_dict = create_label(config.hwdb_test_path[0],config.hwdb_test_path[1])
        data_to_mindrecord_byte_image(test_image_files, test_anno_dict,config.test_dataset_path,prefix="ctpn_test.mindrecord", file_num=1)
    else:
        print("dataset_type should be pretraining, finetune, test")

def data_to_mindrecord_byte_image(image_files, image_anno_dict, dst_dir, prefix="cptn_mlt.mindrecord", file_num=1):
    """Create MindRecord file."""
    mindrecord_path = os.path.join(dst_dir, prefix)
    writer = FileWriter(mindrecord_path, file_num)

    ctpn_json = {
        "image": {"type": "bytes"},
        "annotation": {"type": "int32", "shape": [-1, 5]},
    }
    writer.add_schema(ctpn_json, "ctpn_json")
    for image_name in image_files:
        with open(image_name, 'rb') as f:
            img = f.read()
        annos = np.array(image_anno_dict[image_name], dtype=np.int32)
        print("img name is {}, anno is {}".format(image_name, annos))
        row = {"image": img, "annotation": annos}
        writer.write_raw_data([row])
    writer.commit()

if __name__ == '__main__':
    create_train_dataset("pretraining")
    create_train_dataset("finetune")
    create_train_dataset("test")


text:
import gym
from Agents import QAgent
import numpy as np
from crawler_env import CrawlingRobotEnv

env = CrawlingRobotEnv(render=False)

agent = QAgent(env,gamma=0.9)

all_rewards=0
current_state=env.reset()
total_reward = 0
i = 0
while i < 900000:
	i=i+1
	action = agent.choose_action(current_state)
	next_state,reward,done,info = env.step(action)
	agent.learn(current_state,action,reward,next_state)
	current_state = next_state
	total_reward += reward

	if i % 5000 == 0: # evaluation
	    print "average_reward in last 5000 steps", total_reward / i
	    if (total_reward / i) > 1.3:
	    	break
	    average_reward = 0
	    env.render = True

env = CrawlingRobotEnv(render=True)
current_state=env.reset()
total_reward = 0
agent.eps = 0

i = 0
while True:
	i=i+1
	action = agent.choose_action(current_state)
	next_state,reward,done,info = env.step(action)
	agent.learn(current_state,action,reward,next_state)
	current_state = next_state
	total_reward += reward

	if i % 5000 == 0: # evaluation
	    print "average_reward in last 5000 steps", total_reward / 5000
	    i=0;average_reward = 0
	    env.render = True


text:
#!/usr/bin/env python3

import numpy as np
import rospy
from std_msgs.msg import Float32MultiArray, Bool

from python_op3.framework.modules.walking import Walk

walking_params = [
    "init_x_offset",  # [m]
    "init_y_offset",  # [m]
    "init_z_offset",  # [m]

    "init_roll_offset",  # [radius]
    "init_pitch_offset",  # [radius]
    "init_yaw_offset",  # [radius]

    "hip_pitch_offset",  # [radius]
    "period_time",  # [ms]
    "dsp_ratio",  # [double_stance_time/single_stance_time]

    "step_fb_ratio",  # [m]
    "z_move_amplitude",  # foot_height      [m]
    "y_swap_amplitude",  # swing_right_left [m]

    "z_swap_amplitude",  # swing_top_down   [m]
    "pelvis_offset",  # [radius]
    "arm_swing_gain"  # [m]
]

class AdaptiveWalk(Walk):
    def __init__(self):

        super().__init__()

        self._pub_state = rospy.Publisher("/adaptive_walker/imu_data", Float32MultiArray)
        self._pub_fallen = rospy.Publisher("/adaptive_walker/fallen", Bool)

        self.acc_l = []
        self.gyro_l = []

    def pending_64(self):
        acc_l = []
        gyro_l = []
        while len(acc_l) < 64:
            if self.step_count % 2 == 0:
                acc = self.imu.linear_acceleration
                acc_l.append([acc.x, acc.y, acc.z])
                gyro = self.imu.angular_velocity
                gyro_l.append([gyro.x, gyro.y, gyro.z])
        state = np.array([acc_l, gyro_l])
        return state

    def run_th(self, json_file, ep_i=0):
        self.is_stop = False
        step_params = self._read_ep_params(json_file, ep_i)
        # Set pose
        self._update_walking_params(step_params[0])
        rospy.sleep(5)
        print("json_file:%s epi:%s" % (json_file, ep_i))
        # Enable walking module
        self._pub_module.publish("walking_module")
        self._pub_cmd.publish("start")

        if len(step_params) == 1:
            step_params.extend([[] for _ in range(14)])

        self.step_count = 0
        while not self.is_stop:
            prev_sc = self.step_count
            if prev_sc != self.step_count:
                if len(self.acc_l) < 64:
                    if self.step_count % 2 == 0:
                        acc = self.imu.linear_acceleration
                        self.acc_l.append([acc.x, acc.y, acc.z])
                        gyro = self.imu.angular_velocity
                        self.gyro_l.append([gyro.x, gyro.y, gyro.z])
                    else:
                        rospy.sleep(0.01)
                else:
                    param = self.get_ddpg_predict()
                    self._update_walking_params(param)
            else:
                rospy.sleep(0.001)

        self.stop()

    def get_ddpg_predict(self):
        state = [self.acc_l, self.gyro_l]
        self._pub_state.publish(state)

def run(*args, **kwargs):
    op3 = AdaptiveWalk(*args, **kwargs)
    op3.run()
    # op3._update_walking_params()

if __name__ == '__main__':
    run("up_params.json", epi=0)