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

# # Aprendizado de máquina - II

# ## Resumo
#
# - Na aula anterior... problema de classificação (aprendizado de máquina supervisionado)
# - Nesta aula, resolveremos um problema de _clusterização_ (**aprendizado de máquina não supervisionado**)
# - Aprendizado supervisionado: possível medir erros entre os valores esperados e os previstos
# - Aprendizado não supervisionado: algoritmos categorizam _features_ através de critérios de similaridade, estatísticos ou geométricos.

# ## Breves comentários sobre aprendizagem
#
# - _Aprendizagem_: conceito quantificável por meio de erros.
# - Durante treinamento do algoritmo, dois erros são estimados.

# ### Erro de treinamento
#
# - Mede a discrepância entre as amostras observadas e as amostras apenas do conjunto de treinamento.
# - Erro por amostras (ex. classificação):
# $$e_i = e(x_i,y_i) =
# \begin{cases}
# 1, \ \ \text{se} \ \ h(x_i) = y_i \\
# 0, \ \ \text{caso contrário}
# \end{cases}$$
#
# - Define-se o erro de treinamento por:
# $$E_{tr} = \frac{ \sum_{i=1}^n e_i }{n},$$

# ### Erro de generalização
#
# - Erro médio calculado sobre os dados não vistos ($E_gen$)
# - Em geral, $E_{gen} \geq E_{tr}$ e a _aprendizagem_ é atingida quando $E_{gen}$ é o menor possível.
# - Assim, **a aprendizagem é o resultado de $E_{tr} \to 0$, ou $E_{gen} \approx E_{tr}$**.
# - Um excelente treinamento deve conduzir a um modelo de aprendizagem robusto.

# ### Curvas de aprendizagem, _overfitting_ e _underfitting_
#
# - A aprendizagem é medida por meio de _curvas de aprendizagem_, que possuem aparência assintótica.
# - Modelos com alto desempenho durante a etapa de treinamento que generalizam mal tendem a sofrer um _superajustamento_ (_overfitting_).
# - Modelos com baixo desempenho já durante a etapa de treinamento tendem a sofrer um _subajustamento_ (_underfitting_).
# - Ambos os comportamentos são percebidos nas curvas de aprendizagem.
# - A correção de _overfitting_ ou _underfitting_ pode ser feita com várias técnicas, tais como _regularização_.

# ## Treinamento, validação e generalização
#
# - _treinamento_: etapa em que se verifica o desempenho de um modelo a partir de uma classe de modelos.
# - _validação_: etapa em que se selecionam os melhores hiperparâmetros do modelo tendo por base estimativas do erro de generalização.
# - _generalização (ou teste)_: etapa exclusiva em que se avalia o desempenho final do modelo sobre um conjunto de dados que nunca deve ser usado nas etapas de treinamento ou validação.

# ## Estudo de caso: clusterização de dados aleatórios

# ### O que é clusterização?
#
# Particionamento de dados em subconjuntos disjuntos chamados _clusters_, de forma que os elementos de um cluster sejam altamente similares (alta _similaridade intra-classe_) e os elementos de clusters distintos possuam baixa similaridade com os os demais clusters (baixa _similaridade inter-classe_ ).

# ### Similaridade e distância
#
# - Estabelecer "similaridade" entre instâncias pode ser um pouco difícil.
# - Em problemas de clusterização, usam-se distâncias.

# ### Distância de Minkowski
#
# $$d({\bf x},{\bf y}) = \left( \sum_{i=1}^d | x_i - y_i |^p \right)^{1/p},$$
#
# com ${\bf x}$ e ${\bf y}$ vetores do espaço $\mathbb{R}^n$ e $p \in \mathbb{R}$ é um parâmetro. Dependendo do valor de $p$, a fórmula se reduz a outras distâncias conhecidas:
#
# - Se $p = 2$, temos a _distância Euclidiana_;
# - Se $p = 1$, temos a _distância de Manhattan_ (ou _distância do taxista_);
# - Se $p = \infty$, temos a _distância de Chebyshev_ (ou _distância do máximo_), definida por $\max \{ | x_i - y_i | \}$, isto é, a componente com o maior módulo.

# ### O algoritmo _k-means_
#
# _k-means_ é um algoritmo classificado como de _partição rígida_ que atribui cada elemento no espaço de _features_ a um único cluster. Os passos do algoritmo são:
#
# 1. Inicialize o valor _k_ de clusters desejados;
# 2. Inicialize os centroides dos _k clusters_;
# 3. Atribua cada uma das _n_ amostras de dados ao cluster cujo centroide está mais próximo da amostra;
# 4. Redefina os _k_ centroides iterativamente assumindo que estejam corretos;
# 5. Se nenhuma amostra migrar de um cluster para outro na última iteração, pare. Caso contrário, retorne ao passo 3.

# In[24]:

import numpy as np
import matplotlib.pylab as plt

# In[26]:

# cria dados aleatórios
N = 40

X1 = 1.25*np.random.randn(N,2)
X2 = 5 + 1.5*np.random.randn(N,2)
X3 = [8,3] + 1.2*np.random.randn(N,2)

X = np.concatenate([X1,X2,X3])

X.shape

# Vamos criar rótulos para os 3 grupos e mostrar a clusterização segundo o que nós pré-determinamos e o que, de fato, é passado como entrada para o algoritmo.

# In[28]:

# cria labels
c = np.concatenate([np.full(N,1), np.full(N,2), np.full(N,3)])

plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.scatter(X1[:,0],X1[:,1],color='r')
plt.scatter(X2[:,0],X2[:,1],color='g')
plt.scatter(X3[:,0],X3[:,1],color='b')
plt.title('Dados rotulados')

plt.subplot(1,2,2)
plt.scatter(X[:,0],X[:,1],color='k')
plt.title('Dados "vistos" pelo algoritmo');

# Usa _k-means_ com inicialização de 3 clusters.

# In[30]:

from sklearn import cluster

# no. de clusters
k = 3

# ajusta parâmetros de aprendizagem
cfit = cluster.KMeans(init='random',n_clusters=k)
cfit.fit(X)

# In[32]:

# labels do ajuste
print(cfit.labels_)

# In[35]:

# plota comparativo de labels (original x ajuste)
cl = cfit.labels_ + 1

plt.scatter(X[cl==1,0], X[cl==1,1],color='r')
plt.scatter(X1[:,0],X1[:,1],color='gray',edgecolor='k',alpha=0.4)

plt.scatter(X[cl==2,0], X[cl==2,1],color='g')
plt.scatter(X2[:,0],X2[:,1],color='gray',edgecolor='c',alpha=0.4)

plt.scatter(X[cl==3,0], X[cl==3,1],color='b')
plt.scatter(X3[:,0],X3[:,1],color='gray',edgecolor='m',alpha=0.4);

plt.title('Rótulos originais x rótulos do ajuste');

# Criando pontos novos no espaço que cobre os dados originais.

# In[53]:

# gera "grade" 2D de pontos
# em [-5,15] x [-5,15]
x = np.linspace(-5,15,200)
XX,YY = np.meshgrid(x,x)
sz = XX.shape

# concatena
data = np.c_[XX.flatten(), YY.flatten()]

# predição
Z = cfit.predict(data)

# Verificando partição do espaço determinado pelo algoritmo x dados originais.

# In[56]:

# partições do espaço
plt.imshow(Z.reshape(sz), interpolation='bilinear',
           origin='lower', extent=(-5,15,-5,15),
           alpha=0.3,vmin=0,vmax=k-1)
plt.scatter(X[c==1,0], X[c==1,1],color='r')
plt.scatter(X[c==2,0], X[c==2,1],color='g')
plt.scatter(X[c==3,0], X[c==3,1],color='b')
plt.title('Partições do espaço');

# Então, comparamos o resultado da clusterização em comparação com os dados originais.

# In[62]:

plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
# overlay
plt.imshow(Z.reshape(sz), interpolation='bilinear',
           origin='lower', extent=(-5,15,-5,15),
           alpha=0.3,vmin=0,vmax=k-1)
# original
plt.title('Rótulos originais')
plt.scatter(X[c==1,0], X[c==1,1],color='r')
plt.scatter(X[c==2,0], X[c==2,1],color='g')
plt.scatter(X[c==3,0], X[c==3,1],color='b')
plt.subplot(1,2,2)
# overlay
plt.imshow(Z.reshape(sz), interpolation='bilinear',
           origin='lower', extent=(-5,15,-5,15),
           alpha=0.3,vmin=0,vmax=k-1)
# clusterização
plt.title('Rótulos da clusterização')
plt.scatter(X[cl==1,0], X[cl==1,1],color='r')
plt.scatter(X[cl==2,0], X[cl==2,1],color='g')
plt.scatter(X[cl==3,0], X[cl==3,1],color='b');

# ## Métricas de qualidade para clusterização

# ### Coeficiente de Silhueta
#
# - Valor atribuído a cada ponto de um cluster para medir a "silhueta" de seu cluster
#
# $$s(i) = \frac{b(i) - a(i)}{\max \{a(i),b(i)\}},$$
#
# - $a(i)$ é a distância média intra-cluster da amostra $i$, isto é, a distância média da amostra $i$ a todos as suas "irmãs", percententes ao seu mesmo _cluster_. Este valor mede a dissimilaridade de $i$ em relação às suas amostras "irmãs".
#
# - $b(i)$ é a distância média inter-cluster da amostra $i$, considerando o cluster _A_ ao qual ela pertence e o _cluster_ _B_ mais próximo de _A_. Este valor mede a dissimilaridade de $i$ em relação a todas as amostras em _B_.

# ### Interpretação
#
# O valor de $s$ varia no intervalo [-1,1]
#
# - Quando $s=-1$, a clusterização é considerada incorreta
# - Quando $s=1$, a clusterização está muito densa.
# - Quando $s \approx 0$, há _clusters_ que sobrepõem.

# Vamos calcular $s$ para o nosso exemplo.

# In[66]:

from sklearn import metrics
ss = metrics.silhouette_samples(X,cl,metric='euclidean')
ssm = metrics.silhouette_score(X,cl,metric='euclidean')

# In[69]:

plt.title('Coeficiente de silhueta por amostra')
plt.plot(ss,'+');
plt.plot(np.arange(len(ss)),np.full(np.shape(ss),ssm),':',label='s médio');
plt.legend();

# ### Análise de silhueta
#
# Usada para escolher um valor ótimo para o número de clusters a serem determinados.

# In[71]:

import matplotlib.cm as cm

fig, (ax1,ax2) = plt.subplots(1,2)
fig.set_size_inches(8, 4)

ax1.set_xlim([-0.2,1])
ax1.set_ylim([0, len(X) + (k + 1) * 10])

y_lower = 10
for i in range(1,k+1):
    ssi = ss[cl == i] # i-th cluster s
    ssi.sort()
    sizei = ssi.shape[0]
    y_upper = y_lower + sizei

    color = cm.nipy_spectral(float(i) / k)
    ax1.fill_betweenx(np.arange(y_lower, y_upper),0,
                      ssi,facecolor=color, edgecolor=color,
                      alpha=0.7)
    ax1.text(-0.05, y_lower + 0.5 * sizei, str(i))
    y_lower = y_upper + 10

ax1.set_title(f'Plot de silhueta: {k} clusters')
ax1.set_xlabel("s")
ax1.set_ylabel("Rótulos dos clusters")

ax1.axvline(x=ssm, color="red", linestyle="--")

colors = cm.nipy_spectral(cl.astype(float)/k)
ax2.scatter(X[:, 0],X[:, 1],marker='o', s=30, lw=0, alpha=1.0,
            c=colors, edgecolor='k')
centroides = cfit.cluster_centers_
ax2.scatter(centroides[:, 0], centroides[:, 1], marker='o',
            c="white", alpha=1, s=200, edgecolor='k')

for i, c in enumerate(centroides):

    ax2.scatter(c[0], c[1], marker=f'$%d$' % (i+1), alpha=1,
                s=50, edgecolor='k')

ax2.set_title("Clusters")
ax2.set_xlabel("Feature #1")
ax2.set_ylabel("Feature #2")
plt.suptitle(("Análise de silhueta por k-means: %d clusters" % k),
                 fontsize=12, fontweight='bold');

# ## Conclusão
#
# A análise de silhueta mostra que 3 clusters é um número adequado porque em todos eles há amostras com valor de $s$ acima da média. A espessura das silhuetas também ajuda a ter uma ideia do "tamanho" dos _clusters_.

max_stars_repo_name: gcpeixoto/ICD

max_issues_repo_name: gcpeixoto/ICD

max_forks_repo_name: gcpeixoto/ICD


text:
import unittest

import numpy as np
import numpy.testing as npt
import wisdem.drivetrainse.generator_models as gm

class TestGenerators(unittest.TestCase):
    def setUp(self):
        self.inputs = {}
        self.outputs = {}
        self.discrete_inputs = {}
        self.discrete_outputs = {}

        self.inputs["machine_rating"] = 5e6
        self.inputs["rated_torque"] = 4.143289e6

        self.inputs["rho_Fe"] = 7700.0
        self.inputs["rho_Fes"] = 7850.0
        self.inputs["rho_Copper"] = 8900.0
        self.inputs["rho_PM"] = 7450.0

        self.inputs["B_r"] = 1.2
        self.inputs["E"] = 2e11
        self.inputs["G"] = 79.3e9
        self.inputs["P_Fe0e"] = 1.0
        self.inputs["P_Fe0h"] = 4.0
        self.inputs["S_N"] = -0.002
        self.inputs["alpha_p"] = 0.5 * np.pi * 0.7
        self.inputs["b_r_tau_r"] = 0.45
        self.inputs["b_ro"] = 0.004
        self.inputs["b_s_tau_s"] = 0.45
        self.inputs["b_so"] = 0.004
        self.inputs["cofi"] = 0.85
        self.inputs["freq"] = 60
        self.inputs["h_i"] = 0.001
        self.inputs["h_sy0"] = 0.0
        self.inputs["h_w"] = 0.005
        self.inputs["k_fes"] = 0.9
        self.inputs["k_fillr"] = 0.7
        self.inputs["k_fills"] = 0.65
        self.inputs["k_s"] = 0.2
        self.discrete_inputs["m"] = 3
        self.inputs["mu_0"] = np.pi * 4e-7
        self.inputs["mu_r"] = 1.06
        self.inputs["p"] = 3.0
        self.inputs["phi"] = np.deg2rad(90)
        self.discrete_inputs["q1"] = 6
        self.discrete_inputs["q2"] = 4
        self.inputs["ratio_mw2pp"] = 0.7
        self.inputs["resist_Cu"] = 1.8e-8 * 1.4
        self.inputs["sigma"] = 40e3
        self.inputs["v"] = 0.3
        self.inputs["y_tau_p"] = 1.0
        self.inputs["y_tau_pr"] = 10.0 / 12

    def testPMSG_Outer(self):

        myobj = gm.PMSG_Outer(n_pc=20)

        self.inputs["machine_rating"] = 10.321e6
        self.inputs["rated_torque"] = 10.25e6  # rev 1 9.94718e6
        self.inputs["P_mech"] = 10.71947704e6  # rev 1 9.94718e6
        self.inputs["shaft_rpm"] = np.linspace(2, 10, 20)
        self.inputs["rad_ag"] = 4.0  # rev 1  4.92
        self.inputs["len_s"] = 1.7  # rev 2.3
        self.inputs["h_s"] = 0.7  # rev 1 0.3
        self.inputs["p"] = 70  # 100.0    # rev 1 160
        self.inputs["h_m"] = 0.005  # rev 1 0.034
        self.inputs["h_ys"] = 0.04  # rev 1 0.045
        self.inputs["h_yr"] = 0.06  # rev 1 0.045
        self.inputs["b"] = 2.0
        self.inputs["c"] = 5.0
        self.inputs["B_tmax"] = 1.9
        self.inputs["E_p"] = 3300 / np.sqrt(3)
        self.inputs["D_nose"] = 2 * 1.1  # Nose outer radius
        self.inputs["D_shaft"] = 2 * 1.34  # Shaft outer radius =(2+0.25*2+0.3*2)*0.5
        self.inputs["t_r"] = 0.05  # Rotor disc thickness
        self.inputs["h_sr"] = 0.04  # Rotor cylinder thickness
        self.inputs["t_s"] = 0.053  # Stator disc thickness
        self.inputs["h_ss"] = 0.04  # Stator cylinder thickness
        self.inputs["y_sh"] = 0.0005 * 0  # Shaft deflection
        self.inputs["theta_sh"] = 0.00026 * 0  # Slope at shaft end
        self.inputs["y_bd"] = 0.0005 * 0  # deflection at bedplate
        self.inputs["theta_bd"] = 0.00026 * 0  # Slope at bedplate end
        self.inputs["u_allow_pcent"] = 8.5  # % radial deflection
        self.inputs["y_allow_pcent"] = 1.0  # % axial deflection
        self.inputs["z_allow_deg"] = 0.05  # torsional twist
        self.inputs["sigma"] = 60.0e3  # Shear stress
        self.inputs["B_r"] = 1.279
        self.inputs["ratio_mw2pp"] = 0.8
        self.inputs["h_0"] = 5e-3
        self.inputs["h_w"] = 4e-3
        self.inputs["k_fes"] = 0.8

        myobj.compute(self.inputs, self.outputs, self.discrete_inputs, self.discrete_outputs)
        self.assertAlmostEqual(self.outputs["R_out"], 4.105)
        self.assertAlmostEqual(self.outputs["K_rad"], 0.2125)
        self.assertAlmostEqual(self.outputs["Slot_aspect_ratio"], 6.24779759)
        self.assertAlmostEqual(self.outputs["tau_p"], 0.17951958)
        self.assertAlmostEqual(self.outputs["tau_s"], 0.14930045)
        self.assertAlmostEqual(self.outputs["b_s"], 0.11203948)
        self.assertAlmostEqual(self.outputs["b_t"], 0.03652435)
        self.assertAlmostEqual(self.outputs["h_t"], 0.709)
        self.assertAlmostEqual(self.outputs["b_m"], 0.14361566)
        self.assertAlmostEqual(self.outputs["B_g"], 0.56284833)
        self.assertAlmostEqual(self.outputs["B_symax"], 0.66753916)
        self.assertAlmostEqual(self.outputs["B_rymax"], 0.44502611)
        self.assertAlmostEqual(self.outputs["B_pm1"], 0.46480944)
        self.assertAlmostEqual(self.outputs["B_smax"], 0.23570673147616897)
        self.assertAlmostEqual(self.outputs["f"][-1], 11.66666667)
        self.assertAlmostEqual(self.outputs["I_s"][-1], 1721.4993685170261)
        self.assertAlmostEqual(self.outputs["R_s"][-1], 0.0862357)
        self.assertAlmostEqual(self.outputs["L_s"], 0.01000474)
        self.assertAlmostEqual(self.outputs["S"], 168.0)
        self.assertAlmostEqual(self.outputs["N_s"], 361.0)
        self.assertAlmostEqual(self.outputs["A_Cuscalc"], 389.46615218)
        self.assertAlmostEqual(self.outputs["J_actual"][-1], 3.125519048273891)
        self.assertAlmostEqual(self.outputs["A_1"][-1], 148362.95007673657)
        self.assertAlmostEqual(len(self.outputs["eandm_efficiency"]), 20)
        self.assertAlmostEqual(self.outputs["eandm_efficiency"][-1], 0.9543687904168252)
        self.assertAlmostEqual(self.outputs["Iron"], 89073.14254723)
        self.assertAlmostEqual(self.outputs["mass_PM"], 1274.81620149)
        self.assertAlmostEqual(self.outputs["Copper"], 13859.17179278)
        self.assertAlmostEqual(self.outputs["twist_r"], 0.00032341)
        self.assertAlmostEqual(self.outputs["twist_s"], 5.8057978e-05)
        self.assertAlmostEqual(self.outputs["Structural_mass"], 62323.08483264)
        self.assertAlmostEqual(self.outputs["generator_mass"], 166530.21537414)

    def testPMSG_Arms(self):

        myobj = gm.PMSG_Arms(n_pc=20)

        self.inputs["shaft_rpm"] = np.linspace(5, 12.1, 20)
        self.inputs["sigma"] = 40e3
        self.inputs["rad_ag"] = 3.26
        self.inputs["len_s"] = 1.60
        self.inputs["h_s"] = 0.070
        self.inputs["tau_p"] = 0.080
        self.inputs["h_m"] = 0.009
        self.inputs["h_ys"] = 0.075
        self.inputs["h_yr"] = 0.075
        self.inputs["n_s"] = 5.0
        self.inputs["b_st"] = 0.480
        self.inputs["n_r"] = 5.0
        self.inputs["b_arm"] = 0.530
        self.inputs["d_r"] = 0.700
        self.inputs["d_s"] = 0.350
        self.inputs["t_d"] = 0.105
        self.inputs["t_wr"] = 0.06
        self.inputs["t_ws"] = 0.06
        self.inputs["D_shaft"] = 2 * 0.43
        self.discrete_inputs["q1"] = 1

        myobj.compute(self.inputs, self.outputs, self.discrete_inputs, self.discrete_outputs)
        self.assertAlmostEqual(self.outputs["B_symax"], 0.31869052515138424)
        self.assertAlmostEqual(self.outputs["B_tmax"], 1.391188727217009)
        self.assertAlmostEqual(self.outputs["B_smax"], 0.0995499652108797)
        self.assertAlmostEqual(self.outputs["B_rymax"], 0.28682147263624586)
        self.assertAlmostEqual(self.outputs["B_pm1"], 0.674462389329781)
        self.assertAlmostEqual(self.outputs["B_g"], 0.7651537999693548)
        self.assertAlmostEqual(self.outputs["N_s"], 256.0)
        self.assertAlmostEqual(self.outputs["b_s"], 0.012001865684417254)
        self.assertAlmostEqual(self.outputs["b_t"], 0.01466894694762109)
        self.assertAlmostEqual(self.outputs["A_Cuscalc"], 253.53941258331452)
        self.assertAlmostEqual(self.outputs["b_m"], 0.055999999999999994)
        self.assertAlmostEqual(self.outputs["E_p"][-1], 2013.9466358469522)
        self.assertAlmostEqual(self.outputs["f"][-1], 25.813333333333333)
        self.assertAlmostEqual(self.outputs["I_s"][-1], 846.2065921065657)
        self.assertAlmostEqual(self.outputs["R_s"], 0.09770712863768102)
        self.assertAlmostEqual(self.outputs["L_s"], 0.011401106062847888)
        self.assertAlmostEqual(self.outputs["A_1"][-1], 63455.62872651724)
        self.assertAlmostEqual(self.outputs["J_s"][-1], 3.337574160500578)
        self.assertAlmostEqual(self.outputs["Losses"][-1], 351360.55854497873)
        self.assertAlmostEqual(self.outputs["K_rad"], 0.245398773006135)
        self.assertAlmostEqual(len(self.outputs["eandm_efficiency"]), 20)
        self.assertAlmostEqual(self.outputs["eandm_efficiency"][-1], 0.9343418267745142)
        self.assertAlmostEqual(self.outputs["S"], 768.0)
        self.assertAlmostEqual(self.outputs["Slot_aspect_ratio"], 5.832426544390113)
        self.assertAlmostEqual(self.outputs["Copper"], 6654.6915566955695)
        self.assertAlmostEqual(self.outputs["Iron"], 52673.28069829149)
        self.assertAlmostEqual(self.outputs["TC1"], 16.4856231252069)
        self.assertAlmostEqual(self.outputs["R_out"], 3.3680954773869343)
        self.assertAlmostEqual(self.outputs["Structural_mass"], 33718.05538799999)
        self.assertAlmostEqual(self.outputs["generator_mass"], 94729.99806753898)
        self.assertAlmostEqual(self.outputs["mass_PM"], 1683.970424551947)

    def testPMSG_disc(self):

        myobj = gm.PMSG_Disc(n_pc=20)

        self.inputs["shaft_rpm"] = np.linspace(5, 12.1, 20)
        self.inputs["sigma"] = 40e3
        self.inputs["rad_ag"] = 3.49
        self.inputs["len_s"] = 1.5
        self.inputs["h_s"] = 0.06
        self.inputs["tau_p"] = 0.07
        self.inputs["h_m"] = 0.0105
        self.inputs["h_ys"] = 0.085
        self.inputs["h_yr"] = 0.055
        self.inputs["n_s"] = 5.0
        self.inputs["b_st"] = 0.460
        self.inputs["t_d"] = 0.105
        self.inputs["d_s"] = 0.350
        self.inputs["t_ws"] = 0.150
        self.inputs["D_shaft"] = 2 * 0.43
        self.discrete_inputs["q1"] = 1

        myobj.compute(self.inputs, self.outputs, self.discrete_inputs, self.discrete_outputs)
        self.assertAlmostEqual(self.outputs["B_symax"], 0.2552951)
        self.assertAlmostEqual(self.outputs["B_smax"], 0.0687816)
        self.assertAlmostEqual(self.outputs["B_tmax"], 1.4426429)
        self.assertAlmostEqual(self.outputs["B_rymax"], 0.3550922)
        self.assertAlmostEqual(self.outputs["B_pm1"], 0.6994079)
        self.assertAlmostEqual(self.outputs["B_g"], 0.7934536)
        self.assertAlmostEqual(self.outputs["N_s"], 314.0000000)
        self.assertAlmostEqual(self.outputs["b_s"], 0.0104753)
        self.assertAlmostEqual(self.outputs["b_t"], 0.0128032)
        self.assertAlmostEqual(self.outputs["A_Cuscalc"], 187.2461758)
        self.assertAlmostEqual(self.outputs["b_m"], 0.0490000)
        self.assertAlmostEqual(self.outputs["E_p"][-1], 2555.3443951)
        self.assertAlmostEqual(self.outputs["f"][-1], 31.6616667)
        self.assertAlmostEqual(self.outputs["I_s"][-1], 657.7809841)
        self.assertAlmostEqual(self.outputs["R_s"], 0.1504414)
        self.assertAlmostEqual(self.outputs["L_s"], 0.0120779)
        self.assertAlmostEqual(self.outputs["A_1"][-1], 56514.1132206)
        self.assertAlmostEqual(self.outputs["J_s"][-1], 3.5129208)
        self.assertAlmostEqual(self.outputs["Losses"][-1], 338164.5549178)
        self.assertAlmostEqual(self.outputs["K_rad"], 0.2148997)
        self.assertAlmostEqual(len(self.outputs["eandm_efficiency"]), 20)
        self.assertAlmostEqual(self.outputs["eandm_efficiency"][-1], 0.936651530)
        self.assertAlmostEqual(self.outputs["S"], 942.0000000)
        self.assertAlmostEqual(self.outputs["Slot_aspect_ratio"], 5.7277538)
        self.assertAlmostEqual(self.outputs["Copper"], 5588.6107806)
        self.assertAlmostEqual(self.outputs["Iron"], 48400.5614350)
        self.assertAlmostEqual(self.outputs["TC1"], 16.4856231)
        self.assertAlmostEqual(self.outputs["R_out"], 3.6073317)
        self.assertAlmostEqual(self.outputs["Structural_mass"], 67726.0753421)
        self.assertAlmostEqual(self.outputs["generator_mass"], 123674.5978407)
        self.assertAlmostEqual(self.outputs["mass_PM"], 1959.3502831)

    def testEESG(self):

        myobj = gm.EESG(n_pc=20)

        self.inputs["shaft_rpm"] = np.linspace(5, 12.1, 20)
        self.inputs["sigma"] = 48.373e3
        self.inputs["rad_ag"] = 3.2
        self.inputs["len_s"] = 1.4
        self.inputs["h_s"] = 0.060
        self.inputs["tau_p"] = 0.170
        self.inputs["I_f"] = 69
        self.inputs["N_f"] = 100
        self.inputs["h_ys"] = 0.130
        self.inputs["h_yr"] = 0.120
        self.inputs["n_s"] = 5
        self.inputs["b_st"] = 0.470
        self.inputs["n_r"] = 5
        self.inputs["b_arm"] = 0.480
        self.inputs["d_r"] = 0.510
        self.inputs["d_s"] = 0.400
        self.inputs["t_wr"] = 0.140
        self.inputs["t_ws"] = 0.070
        self.inputs["D_shaft"] = 2 * 0.43
        self.discrete_inputs["q1"] = 2

        myobj.compute(self.inputs, self.outputs, self.discrete_inputs, self.discrete_outputs)
        self.assertAlmostEqual(self.outputs["B_symax"], 0.45940641371625035)
        self.assertAlmostEqual(self.outputs["B_tmax"], 1.8877594585933197)
        self.assertAlmostEqual(self.outputs["B_rymax"], 0.6181179999135469)
        self.assertAlmostEqual(self.outputs["B_g"], 1.1036753874387688)
        self.assertAlmostEqual(self.outputs["N_s"], 236.0)
        self.assertAlmostEqual(self.outputs["b_s"], 0.012779359946805938)
        self.assertAlmostEqual(self.outputs["b_t"], 0.015619217712762814)
        self.assertAlmostEqual(self.outputs["A_Cuscalc"], 228.43105904915618)
        self.assertAlmostEqual(self.outputs["f"][-1], 11.898333333333333)
        self.assertAlmostEqual(self.outputs["I_s"][-1], 858.2107484477335)
        self.assertAlmostEqual(self.outputs["R_s"], 0.0906017600502669)
        self.assertAlmostEqual(self.outputs["A_1"][-1], 60440.40365229799)
        self.assertAlmostEqual(self.outputs["J_s"][-1], 3.75697924800609)
        self.assertAlmostEqual(self.outputs["Losses"][-1], 444556.7203870894)
        self.assertAlmostEqual(self.outputs["K_rad"], 0.21874999999999997)
        self.assertAlmostEqual(len(self.outputs["eandm_efficiency"]), 20)
        self.assertAlmostEqual(self.outputs["eandm_efficiency"][-1], 0.918348408655116)
        self.assertAlmostEqual(self.outputs["S"], 708.0)
        self.assertAlmostEqual(self.outputs["Slot_aspect_ratio"], 4.695070821210912)
        self.assertAlmostEqual(self.outputs["Copper"], 16059.414075335235)
        self.assertAlmostEqual(self.outputs["Iron"], 72360.88485924648)
        self.assertAlmostEqual(self.outputs["TC1"], 13.632086598066605)
        self.assertAlmostEqual(self.outputs["R_out"], 3.219748743718593)
        self.assertAlmostEqual(self.outputs["Structural_mass"], 42403.44234)
        self.assertAlmostEqual(self.outputs["generator_mass"], 130823.74127458173)

    def testSCIG(self):

        myobj = gm.SCIG(n_pc=20)

        self.inputs["shaft_rpm"] = 100 * np.linspace(2, 12, 20)
        self.inputs["cofi"] = 0.9
        self.inputs["y_tau_p"] = 12.0 / 15.0
        self.inputs["sigma"] = 21.5e3
        self.inputs["rad_ag"] = 0.55
        self.inputs["len_s"] = 1.30
        self.inputs["h_s"] = 0.090
        self.inputs["I_0"] = 140
        self.inputs["B_symax"] = 1.4
        self.discrete_inputs["q1"] = 6
        self.inputs["h_0"] = 0.05

        myobj.compute(self.inputs, self.outputs, self.discrete_inputs, self.discrete_outputs)
        self.assertAlmostEqual(self.outputs["B_tsmax"], 1.684969767820678)
        self.assertAlmostEqual(self.outputs["B_trmax"], 1.8541950805798264)
        self.assertAlmostEqual(self.outputs["B_rymax"], 1.4)
        self.assertAlmostEqual(self.outputs["B_g"], 0.9267333723013729)
        self.assertAlmostEqual(self.outputs["B_g1"], 0.8617526608763629)
        self.assertAlmostEqual(self.outputs["h_ys"], 0.12135794161089408)
        self.assertAlmostEqual(self.outputs["b_s"], 0.014398966328953218)
        self.assertAlmostEqual(self.outputs["b_t"], 0.017598736624276155)
        self.assertAlmostEqual(self.outputs["D_ratio"], 1.393378075656171)
        self.assertAlmostEqual(self.outputs["D_ratio_UL"], 1.4)
        self.assertAlmostEqual(self.outputs["D_ratio_LL"], 1.37)
        self.assertAlmostEqual(self.outputs["A_Cuscalc"], 244.78242759220464)
        self.assertAlmostEqual(self.outputs["Slot_aspect_ratio1"], 6.250448674154435)
        self.assertAlmostEqual(self.outputs["h_yr"], 0.12135794161089408)
        self.assertAlmostEqual(self.outputs["tau_p"], 0.5759586531581288)
        self.assertAlmostEqual(self.outputs["b_r"], 0.021513941784534422)
        self.assertAlmostEqual(self.outputs["b_trmin"], 0.023894990015160978)
        self.assertAlmostEqual(self.outputs["b_tr"], 0.026294817736653178)
        self.assertAlmostEqual(self.outputs["rad_r"], 0.547848028863988)
        self.assertAlmostEqual(self.outputs["A_bar"], 0.000968127380304049)
        self.assertAlmostEqual(self.outputs["Slot_aspect_ratio2"], 2.3240743375043973)
        self.assertAlmostEqual(self.outputs["E_p"][-1], 3584.6363848734995)
        self.assertAlmostEqual(self.outputs["I_s"][-1], 383.73299436233305)
        self.assertAlmostEqual(self.outputs["A_1"][-1], 23985.033858413553)
        self.assertAlmostEqual(self.outputs["J_s"][-1], 1.5676492717917363)
        self.assertAlmostEqual(self.outputs["J_r"][-1], 0.3690452555761286)
        self.assertAlmostEqual(self.outputs["R_s"], 0.00820627176284264)
        self.assertAlmostEqual(self.outputs["R_R"], 0.026060741300117953)
        self.assertAlmostEqual(self.outputs["L_s"], 0.001078616172170861)
        self.assertAlmostEqual(self.outputs["L_sm"], 0.06791816664758084)
        self.assertAlmostEqual(self.outputs["generator_mass"], 40729.31598698678)
        self.assertAlmostEqual(self.outputs["K_rad"], 1.1818181818181817)
        self.assertAlmostEqual(self.outputs["Losses"][-1], 75068.46569051298)
        self.assertAlmostEqual(len(self.outputs["eandm_efficiency"]), 20)
        self.assertAlmostEqual(self.outputs["eandm_efficiency"][-1], 0.9849562794756211)
        self.assertAlmostEqual(self.outputs["Copper"], 1360.253987369221)
        self.assertAlmostEqual(self.outputs["Iron"], 11520.904698026085)
        self.assertAlmostEqual(self.outputs["R_out"], 0.7663579416108941)
        self.assertAlmostEqual(self.outputs["Structural_mass"], 27848.15730159148)
        self.assertAlmostEqual(self.outputs["TC1"], 0.29453832454167966)

    def testDFIG(self):

        myobj = gm.DFIG(n_pc=20)

        self.inputs["shaft_rpm"] = 100 * np.linspace(2, 12, 20)
        self.inputs["cofi"] = 0.9
        self.inputs["y_tau_p"] = 12.0 / 15.0
        self.inputs["sigma"] = 21.5e3
        self.inputs["rad_ag"] = 0.61
        self.inputs["len_s"] = 0.49
        self.inputs["h_s"] = 0.08
        self.inputs["I_0"] = 40.0
        self.inputs["B_symax"] = 1.3
        self.inputs["S_Nmax"] = -0.2
        self.inputs["k_fillr"] = 0.55
        self.inputs["h_0"] = 0.01
        self.inputs["k_fillr"] = 0.55
        self.inputs["k_fills"] = 0.65
        self.discrete_inputs["q1"] = 5
        self.discrete_inputs["q2"] = 4

        myobj.compute(self.inputs, self.outputs, self.discrete_inputs, self.discrete_outputs)
        self.assertAlmostEqual(self.outputs["B_g"], 1.1073196940043413)
        self.assertAlmostEqual(self.outputs["B_g1"], 1.0207571181956336)
        self.assertAlmostEqual(self.outputs["B_rymax"], 1.3)
        self.assertAlmostEqual(self.outputs["B_tsmax"], 2.0133085345533477)
        self.assertAlmostEqual(self.outputs["B_trmax"], 2.046984459493319)
        self.assertAlmostEqual(self.outputs["h_ys"], 0.1731961572673457)
        self.assertAlmostEqual(self.outputs["b_s"], 0.019163715186897738)
        self.assertAlmostEqual(self.outputs["b_t"], 0.023422318561763904)
        self.assertAlmostEqual(self.outputs["D_ratio"], 1.4232723889628618)
        self.assertAlmostEqual(self.outputs["A_Cuscalc"], 287.4557278034661)
        self.assertAlmostEqual(self.outputs["Slot_aspect_ratio1"], 4.1745558843775825)
        self.assertAlmostEqual(self.outputs["h_yr"], 0.1731961572673457)
        self.assertAlmostEqual(self.outputs["tau_p"], 0.6387905062299246)
        self.assertAlmostEqual(self.outputs["N_r"], 148.0)
        self.assertAlmostEqual(self.outputs["b_r"], 0.023870136274726762)
        self.assertAlmostEqual(self.outputs["b_trmin"], 0.02869464545814538)
        self.assertAlmostEqual(self.outputs["b_tr"], 0.02917461100244382)
        self.assertAlmostEqual(self.outputs["A_Curcalc"], 5.322395250445834)
        self.assertAlmostEqual(self.outputs["Slot_aspect_ratio2"], 0.41893351109971694)
        self.assertAlmostEqual(self.outputs["E_p"][-1], 1480.0153123457287)
        self.assertAlmostEqual(self.outputs["I_s"][-1], 959.242520719389)
        self.assertAlmostEqual(self.outputs["A_1"][-1], 45049.62948091099)
        self.assertAlmostEqual(self.outputs["J_s"][-1], 3.3370095911785915)
        self.assertAlmostEqual(self.outputs["J_r"][-1], 176.31699359563828)
        self.assertAlmostEqual(self.outputs["R_s"], 0.0037830355111420555)
        self.assertAlmostEqual(self.outputs["R_R"], 0.03168834286109277)
        self.assertAlmostEqual(self.outputs["L_r"], 0.00016404053494963115)
        self.assertAlmostEqual(self.outputs["L_s"], 0.0001855520925628366)
        self.assertAlmostEqual(self.outputs["L_sm"], 0.019568172828126774)
        self.assertAlmostEqual(self.outputs["generator_mass"], 19508.405570203206)
        self.assertAlmostEqual(self.outputs["K_rad"], 0.4016393442622951)
        self.assertAlmostEqual(self.outputs["Losses"][-1], 143901.77088462992)
        self.assertAlmostEqual(len(self.outputs["eandm_efficiency"]), 20)
        self.assertAlmostEqual(self.outputs["eandm_efficiency"][-1], 0.9654635749876888)
        self.assertAlmostEqual(self.outputs["Copper"], 354.97950026786026)
        self.assertAlmostEqual(self.outputs["Iron"], 6077.944572132577)
        self.assertAlmostEqual(self.outputs["Structural_mass"], 13075.48149780277)
        self.assertAlmostEqual(self.outputs["R_out"], 0.8681961572673457)
        self.assertAlmostEqual(self.outputs["TC1"], 0.24790308982258036)
        self.assertAlmostEqual(self.outputs["Current_ratio"][-1], 0.17760183748148742)

def suite():
    suite = unittest.TestSuite()
    suite.addTest(unittest.makeSuite(TestGenerators))
    return suite

if __name__ == "__main__":
    result = unittest.TextTestRunner().run(suite())

    if result.wasSuccessful():
        exit(0)
    else:
        exit(1)

ext: py

max_stars_repo_name: johnjasa/WISDEM

max_issues_repo_name: johnjasa/WISDEM

max_forks_repo_name: johnjasa/WISDEM

reason: import numpy


text:
# -*- tab-width:4;indent-tabs-mode:nil;show-trailing-whitespace:t;rm-trailing-spaces:t -*-
# vi: set ts=2 noet:

import time
import os
import shutil
import joblib
import numpy as np
import pyarrow.parquet
import pyarrow as pa
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import (PCA, IncrementalPCA)
import umap
import pandas as pd
import datashader
from datashader.transfer_functions import set_background
from colorcet import fire

def fit_embedding(
        dataset,
        embed_dir,
        standardize_features=True,
        pca_n_components=None,
        pca_batch_size=1000,
        umap_n_components=2,
        umap_init='random',
        umap_n_neighbors=100,
        umap_min_dist=0.0,
        umap_a=None,
        umap_b=None,
        umap_negative_sample_rate=5,
        umap_metric='euclidean',
        umap_n_epochs=None,
        low_memory=False,
        save_transform=True,
        seed=None,
        verbose=True):
    """
    train_set: a feature matrix e.g. of (N, F) dimensions, used to define the embedding
    After some experimentation with chemical features of 1024 of 50k-.5M compounds,
    a reasonable embedding first reduces by PCA to 20 dimensions and then UMAP to 2 dimensions.
    UMAP parameters of 100 neighbors and min_dist of 0.0 seem to work well too. init='random'
    can help UMAP from getting stuck.
    return:
        saves embedding data to
            ../intermediate_data/embeddings/tag/embedding_info.tsv
            ../intermediate_data/embeddings/tag/pca_reducer.joblib
            ../intermediate_data/embeddings/tag/umap_reducer.joblib
    """
    if not os.path.exists(embed_dir):
        os.mkdir(embed_dir)
    else:
        print("WARNING: embed_dir already exists: {}".format(embed_dir))

    random_state = np.random.RandomState(seed=seed)

    begin_time = time.time()

    if standardize_features:
        if verbose:
            print("Standardizing dataset so each feature has zero-mean and unit variance.")
        standardizer = StandardScaler(copy=False)
        standardizer.fit(dataset)
        dataset = standardizer.transform(dataset)

    if pca_n_components is None:
        pca_n_components = min(dataset.shape[0], dataset.shape[1])
        if verbose:
            print("Setting PCA n_componets to full rank of dataset: {}".format(pca_n_components))
    elif pca_n_components > dataset.shape[1]:
        pca_n_components = dataset.shape[1]
        if verbose:
            print(f"Requesting PCA n_components of {pca_n_components} but the dataset only has {dataset.shape[1]} columns, setting pca_n_components to {dataset.shape[1]}")
    elif pca_n_components > dataset.shape[0]:
        pca_n_components = dataset.shape[0]
        if verbose:
            print(f"Requesting PCA n_components of {pca_n_components} but the dataset only has {dataset.shape[0]} rows, setting pca_n_components to {dataset.shape[0]}")

    if pca_batch_size < pca_n_components:
        if verbose:
            print(f"Adjusting PCA batch size to be at least the number of components '{pca_n_components}'")
        pca_batch_size = pca_n_components

    if verbose:
        print("Reducing the dimension by PCA from {} to {} dimensions".format(dataset.shape[1], pca_n_components))
    pca_reducer = IncrementalPCA(
        n_components=pca_n_components,
        batch_size=pca_batch_size,
        copy=False)
    pca_reducer.fit(dataset)
    pca_embedding = pca_reducer.transform(dataset)

    if verbose:
        print("Reducing the dimension by UMAP to {} dimensions".format(umap_n_components))
    umap_reducer = umap.UMAP(
        n_components=umap_n_components,
        metric=umap_metric,
        n_neighbors=umap_n_neighbors,
        min_dist=umap_min_dist,
        a=umap_a,
        b=umap_b,
        negative_sample_rate=umap_negative_sample_rate,
        init=umap_init,
        low_memory=low_memory,
        random_state=random_state,
        n_epochs=umap_n_epochs,
        verbose=True)
    umap_embedding = umap_reducer.fit_transform(pca_embedding)
    umap_embedding = pd.DataFrame(
        data=umap_embedding,
        columns=["UMAP_" + str(i+1) for i in range(umap_n_components)])
    end_time = time.time()

    if verbose:
        print("created embedding {0} runtime: {1:.2f}s".format(embed_dir, end_time-begin_time))
        print("saving embedding to {}".format(embed_dir))
    with open("{}/model_info.tsv".format(embed_dir), 'w') as f:
        f.write("key\tvalue\n")
        f.write("seed\t{}\n".format(seed))
        f.write("input_dim\t{}\n".format(dataset.shape))
        f.write("standardize_features\t{}\n".format(standardize_features))
        f.write("pca_n_component\t{}\n".format(pca_n_components))
        f.write("umap_n_component\t{}\n".format(umap_n_components))
        f.write("umap_metric\t{}\n".format(umap_metric))
        f.write("umap_n_neighbors\t{}\n".format(umap_n_neighbors))
        f.write("umap_min_dist\t{}\n".format(umap_min_dist))
        f.write("umap_a\t{}\n".format(umap_reducer._a))
        f.write("umap_b\t{}\n".format(umap_reducer._b))
        f.write("umap_negative_sample_rate\t{}\n".format(umap_negative_sample_rate))
        f.write("umap_low_memory\t{}\n".format(low_memory))
        f.write("umap_init\t{}\n".format(umap_init))
        f.write("umap_n_epochs\t{}\n".format(umap_n_epochs))

    pa.parquet.write_table(
        table=pa.Table.from_pandas(umap_embedding),
        where="{}/umap_embedding.parquet".format(embed_dir))

    if save_transform:
        if verbose:
            print("Saving transform to {}.".format(embed_dir))
        if pca_n_components is not None:
            joblib.dump(
                value=pca_reducer,
                filename="{}/pca_reducer.joblib".format(embed_dir))
        if standardize_features:
            joblib.dump(
                value=standardizer,
                filename="{}/standardizer.joblib".format(embed_dir))
        joblib.dump(
            value=umap_reducer,
            filename="{}/umap_reducer.joblib".format(embed_dir))

    return umap_embedding

def embed(
        dataset,
        embed_dir,
        ref_embed_dir,
        standardize_features=None,
        batch_size=None,
        n_epochs=None,
        verbose=True):
    """
    Given a previously defined embedding, embed a new dataset into it
    """

    if not os.path.exists(embed_dir):
        os.mkdir(embed_dir)
    else:
        print("WARNING: embed_dir already exists: {}".format(embed_dir))

    begin_time = time.time()
    if standardize_features:
        if verbose:
            print("Standardizing dataset so each feature has zero-mean and unit variance.")
        standardizer = StandardScaler(copy=False)
        standardizer.fit(dataset)
        dataset = standardizer.transform(dataset)
    elif standardize_features is None:
        if verbose:
            print("Loading standardizer from embed_dir {} ...".format(ref_embed_dir))
        try:
            standardizer = joblib.load(filename="{}/standardizer.joblib".format(ref_embed_dir))
        except:
            standardizer = None
            if verbose:
                print("Standardizer not found at '{}/standardizer.joblib', assume no-prestandardization is needed".format(ref_embed_dir))
    else:
        standardizer = None

    if verbose:
        print("Loading reference PCA->UMAP reducer from embed_dir {} ...".format(ref_embed_dir))
    try:
        pca_reducer = joblib.load(filename="{}/pca_reducer.joblib".format(ref_embed_dir))
    except:
        if verbose:
            print("PCA Reducer not found at '{}/pca_reducer.joblib', assume no pre-dimensionality reduction by PCA".format(ref_embed_dir))
        pca_reducer = None

    umap_reducer = joblib.load(filename="{}/umap_reducer.joblib".format(ref_embed_dir))
    umap_reducer.n_epochs = n_epochs

    if batch_size is None:
        n_batches = 1
    else:
        n_batches = np.floor(dataset.shape[0]/batch_size)
    umap_embedding = []
    for batch_index, batch in enumerate(np.array_split(dataset, n_batches)):
        if verbose:
            print("Transforming batch {} of {}".format(batch_index+1, n_batches))
        if standardizer:
            batch = standardizer.transform(batch)
        if pca_reducer:
            pca_embedding = pca_reducer.transform(batch)
        else:
            pca_embdding = batch
        umap_embedding.append(umap_reducer.transform(pca_embedding))
    umap_embedding = np.vstack(umap_embedding)
    umap_embedding = pd.DataFrame(
        data=umap_embedding,
        columns=["UMAP_" + str(i+1) for i in range(umap_embedding.shape[1])])
    end_time = time.time()
    if verbose:
        print("embedded data onto {0} runtime: {1:.2f}s".format(embed_dir, end_time-begin_time))

    if verbose:
        print("Copying model info from reference embed dir.")

    try:
        shutil.copyfile(
            src="{}/model_info.tsv".format(ref_embed_dir),
            dst="{}/model_info.tsv".format(embed_dir))
    except:
        print("WARNING: Failed to copy reference model info from '{}/model_info.tsv' to '{}/model_info.tsv'".format(ref_embed_dir, embed_dir))

    with open("{}/model_info.tsv".format(embed_dir), 'a') as f:
        f.write("ref_embed_dir\t{}\n".format(ref_embed_dir))
        f.write("n_epochs\t{}\n".format(n_epochs))

    pa.parquet.write_table(
        table=pa.Table.from_pandas(umap_embedding),
        where="{}/umap_embedding.parquet".format(embed_dir))
    return umap_embedding

max_stars_repo_name: momeara/MPLearn

max_issues_repo_name: momeara/MPLearn

max_forks_repo_name: momeara/MPLearn


text:
from __future__ import print_function

import dynamic_graph.sot_talos_balance.talos.control_manager_conf as control_manager_conf
import dynamic_graph.sot_talos_balance.talos.parameter_server_conf as param_server_conf
import numpy as np
from dynamic_graph.sot_talos_balance.create_entities_utils import ParameterServer, TalosControlManager
from dynamic_graph.sot_talos_balance.utils.sot_utils import Bunch

dt = 0.001
conf = Bunch()
robot_name = 'robot'
N_JOINTS = 38
u_max = 1.

conf.param_server = param_server_conf
conf.control_manager = control_manager_conf

param_server = ParameterServer("param_server")
param_server.init(dt, conf.param_server.urdfFileName, robot_name)
param_server.setJointsUrdfToSot(np.array(conf.param_server.urdftosot))
param_server.setRightFootForceSensorXYZ(np.array(conf.param_server.rightFootSensorXYZ))
param_server.setRightFootSoleXYZ(np.array(conf.param_server.rightFootSoleXYZ))
for key in conf.param_server.mapJointNameToID:
    param_server.setNameToId(key, conf.param_server.mapJointNameToID[key])

cm = TalosControlManager("TalosControlManager")
cm.add_commands()
cm.init(dt, robot_name)
print("***Control manager initialized***")
cm.add_signals()
cm.u_max.value = np.array(N_JOINTS * (u_max, ))
print("u_max is set at {0} for all joints".format(u_max, ))

# Control is Ok
print("*****************")
control_value = 0.5
u = np.array(N_JOINTS * [control_value])
cm.addCtrlMode('input')
cm.setCtrlMode('all', 'input')
cm.add_signals()
cm.ctrl_input.value = u
print("Control input is set at {0} for all joints".format(control_value))
cm.u.recompute(1)
cm.u_safe.recompute(1)
assert (cm.u_safe.value == (control_value, ) * N_JOINTS).all()
print("Safe control = control input")
print("*****************")
# Control is too big
control_value = 2.
u = np.array(N_JOINTS * [control_value])
cm.ctrl_input.value = u
print("Control input is set at {0} for all joints".format(control_value))
cm.u.recompute(2)
cm.u_safe.recompute(2)
assert (cm.u_safe.value == (0., ) * N_JOINTS).all()
print("Safe control = zero")
print("Control set to 0 forever")
print("*****************")

cm2 = TalosControlManager("TalosControlManager")
cm2.add_commands()
cm2.init(dt, robot_name)
print("***New Control manager initialized***")
cm2.add_signals()
cm.u_max.value = np.array(N_JOINTS * (u_max, ))
print("u_max is set at {0} for all joints".format(u_max))

# Control is Ok
print("*****************")
control_value = 0.5
u = np.array(N_JOINTS * [control_value])
cm2.addCtrlMode('input')
cm2.setCtrlMode('all', 'input')
cm2.ctrl_input.value = u
print("Control input is set at {0} for all joints".format(control_value))
cm2.u.recompute(3)
cm2.u_safe.recompute(3)
assert (cm2.u_safe.value == (control_value, ) * N_JOINTS).all()
print("Safe control = control input")
print("*****************")

# Fake emergency signal
emergency = 0
print("EMERGENCY = 0")
cm2.addEmergencyStopSIN('test')
cm2.add_signals()
cm2.emergencyStop_test.value = emergency
cm2.u.recompute(4)
cm2.u_safe.recompute(4)
assert (cm2.u_safe.value == (control_value, ) * N_JOINTS).all()
print("Safe control = control input")
emergency = 1
print('EMERGENCY = 1')
cm2.emergencyStop_test.value = emergency
cm2.u.recompute(5)
cm2.u_safe.recompute(5)
assert (cm2.u_safe.value == (0., ) * N_JOINTS).all()
print("Safe control = zero")
print("Control set to 0 forever")
print("*****************")


text:
# This script combines the intermediate data into processed data
# outputs:
#   node.pickle - shape: (110 * _) - all the nodal information
#   edge:pickle - shape: (12100 * _) - all the edge information

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
import pickle
import copy

# system path
import sys
import os

# util path
utility_path = os.path.join(os.getcwd(),'src/d00_utils/')
sys.path.append(utility_path)
import utilities as util

# data path
# sw: define the path based on the root project directory.
intermediate_data_path = os.path.join(os.getcwd(),'data/02_intermediate/')
processing_data_path = os.path.join(os.getcwd(),'data/03_processed/')

# read files
sa2_adelaide_node = gpd.read_file(intermediate_data_path+'shapefiles/sa2_adelaide.shp')
sa2_adelaide_edge = gpd.read_file(intermediate_data_path+'shapefiles/sa2_adelaide_edge.shp')
# read six files
with open(intermediate_data_path+'sa2_node_with_socio_econ_df.pickle', 'rb') as f:
    sa2_node_with_socio_econ_df = pickle.load(f)

with open(intermediate_data_path+'sa2_node_with_POI_counts_entropy.pickle', 'rb') as f:
    sa2_node_with_POI_counts_entropy = pickle.load(f)

with open(intermediate_data_path+'sa2_node_with_only_transport_attributes.pickle', 'rb') as f:
    sa2_node_with_only_transport_attributes = pickle.load(f)

with open(intermediate_data_path+'sa2_edge_flow.pickle', 'rb') as f:
    sa2_edge_flow = pickle.load(f)

with open(intermediate_data_path+'sa2_edge_travel_time.pickle', 'rb') as f:
    sa2_edge_travel_time = pickle.load(f)

with open(intermediate_data_path+'sa2_edge_with_only_transport_attributes.pickle', 'rb') as f:
    sa2_edge_with_only_transport_attributes = pickle.load(f)

# drop geometry columns
sa2_edge_flow.drop(columns=['geometry'],inplace=True)
sa2_edge_with_only_transport_attributes.drop(columns=['geometry'],inplace=True)
sa2_node_with_POI_counts_entropy.drop(columns=['geometry'],inplace=True)
sa2_node_with_only_transport_attributes.drop(columns=['geometry'],inplace=True)

# only subset of socio-econ variables are needed for now.
sa2_node_with_socio_econ_df = sa2_node_with_socio_econ_df[['sa2_code16',
                                                           'num_jobs_000_persons',
                                                           'median_income_per_job_aud_persons',
                                                           'tot_num_jobs_000',
                                                           'total_pop',
                                                           'avg_age',
                                                           'male_percent',
                                                           'female_percent',
                                                           'bachelor_degree_percent',
                                                           'master_degree_percent',
                                                           'perc_indig_age_0_14',
                                                           'perc_indig_age_15_34',
                                                           'perc_indig_age_35_64',
                                                           'poverty_rate_1',
                                                           'poverty_rate_2',
                                                           'unemployment_rate',
                                                           'median_inc',
                                                           'gini']]

#region 1. merge edge files
# sa2_adelaide_edge, sa2_edge_flow, sa2_edge_travel_time, sa2_edge_with_only_transport_attributes
# disaggregate nodal files: sa2_node_with_POI_counts_entropy, sa2_node_with_socio_econ_df.
edge = copy.copy(sa2_adelaide_edge)

# merge with edge files
edge = edge.merge(sa2_edge_flow, left_on=['O','D'], right_on=['O','D'], how='left')
edge = edge.merge(sa2_edge_travel_time, left_on=['O','D'], right_on=['O','D'], how='left')
edge = edge.merge(sa2_edge_with_only_transport_attributes, left_on=['O','D'], right_on=['O','D'], how='left')

# add nodal files to the edge files
edge = edge.merge(sa2_node_with_POI_counts_entropy, left_on=['O'], right_on=['SA2_MAIN16'], how='left')
edge = edge.merge(sa2_node_with_POI_counts_entropy, left_on=['D'], right_on=['SA2_MAIN16'], suffixes=['_o', '_d'], how='left')
edge = edge.merge(sa2_node_with_socio_econ_df, left_on=['O'], right_on=['sa2_code16'], how='left')
edge = edge.merge(sa2_node_with_socio_econ_df, left_on=['D'], right_on=['sa2_code16'], suffixes=['_o', '_d'], how='left')

# save edge
edge.to_pickle(processing_data_path+'edge_shp.pickle')

#endregion

#region 2. merge nodal files
# sa2_adelaide_node, sa2_node_with_socio_econ_df, sa2_node_with_POI_counts_entropy, sa2_node_with_only_transport_attributes;
node = copy.copy(sa2_adelaide_node)
node = node.merge(sa2_node_with_socio_econ_df, left_on=['SA2_MAIN16'], right_on=['sa2_code16'], how='left')
node = node.merge(sa2_node_with_POI_counts_entropy, left_on=['SA2_MAIN16'], right_on=['SA2_MAIN16'], how='left')
node = node.merge(sa2_node_with_only_transport_attributes, left_on=['SA2_MAIN16'], right_on=['SA2_MAIN16'], how='left')

# aggregate edge info to O and D. sa2_edge_flow
# Then add O and D info to the nodal file
sa2_flow_o = sa2_edge_flow.groupby('O').aggregate(['sum']).reset_index()
col_list = sa2_flow_o.columns
col_new = [col[0]+'_o' for col in col_list]
sa2_flow_o.columns = col_new

sa2_flow_d = sa2_edge_flow.groupby('D').aggregate(['sum']).reset_index()
col_list = sa2_flow_d.columns
col_new = [col[0]+'_d' for col in col_list]
sa2_flow_d.columns = col_new

node = node.merge(sa2_flow_o, left_on=['SA2_MAIN16'], right_on=['O_o'], how='left')
node = node.merge(sa2_flow_d, left_on=['SA2_MAIN16'], right_on=['D_d'], how='left')

# drop two columns
node = node.drop(columns=['D_o', 'O_d'])

# change data types
node[['num_nodes', 'num_1degree', 'num_2degree', 'num_3degree', 'num_4degree', 'num_greater5degree']] = \
    node[['num_nodes', 'num_1degree', 'num_2degree', 'num_3degree', 'num_4degree', 'num_greater5degree']].astype('int64')

# add pop density variable
node['pop_density'] = node['total_pop']/(node.area * 1000)

# add POI densities
node['poi_count_density']=node['poi_count']/(1000*node.area)
node['poi_count_agg_density']=node['poi_count_agg']/(1000*node.area)
node['poi_entropy_density']=node['poi_entropy']/(1000*node.area)
node['poi_entropy_agg_density']=node['poi_entropy_agg']/(1000*node.area)

# add job densities
node['job_density_1']=node['num_jobs_000_persons']/(1000*node.area)
node['job_density_2']=node['tot_num_jobs_000']/(1000*node.area)

# add o and d info together
node['flow_agents'] = node['flow_agents_o'] + node['flow_agents_d']
node['flow_duration'] = node['flow_duration_o'] + node['flow_duration_d']
node['flow_stays'] = node['flow_stays_o'] + node['flow_stays_d']
node['consumption_count_age_source'] = node['consumption_count_age_source_o'] + node['consumption_count_age_source_d']
node['consumption_amount_age_source'] = node['consumption_amount_age_source_o'] + node['consumption_amount_age_source_d']
node['consumption_count_mcc_source'] = node['consumption_count_mcc_source_o'] + node['consumption_count_mcc_source_d']
node['consumption_amount_mcc_source'] = node['consumption_amount_mcc_source_o'] + node['consumption_amount_mcc_source_d']

# save node
node.to_pickle(processing_data_path+'node_shp.pickle')
#endregion

#region 3. processing node and edge shapefiles into three data frames
# read node and edge pickles
with open(processing_data_path+'node_shp.pickle', 'rb') as f:
    node = pickle.load(f)

with open(processing_data_path+'edge_shp.pickle', 'rb') as f:
    edge = pickle.load(f)

# create dataframes to facilitate viewing as dataframes in pycharm
node_df = pd.DataFrame(node.drop(columns=['geometry']))
edge_df = pd.DataFrame(edge.drop(columns=['geometry']))

# edit node_df
print(np.sum(node_df.isna()))
print(np.sum(node_df==0))
print(node_df.describe())
# processing na and zeros for regression on the log space
node_df.dropna(how='any',inplace=True)
print(node_df.shape) # dropped 6 observations. Left with 104 obs.
non_zero_masks = np.logical_and(node_df.bachelor_degree_percent>0.00001, node_df.flow_agents_o>0.00001)
non_zero_masks = np.logical_and(node_df.poi_entropy_agg>0.00001, non_zero_masks)
node_df = node_df.loc[non_zero_masks, :]
print(node_df.shape) # dropped 2 observations. Left with 102 obs.

# lift all road attributes by one for the log transformation.
road_attributes = ['class_ART',
                   'class_BUS', 'class_COLL', 'class_HWY', 'class_LOCL',
                   'class_SUBA', 'class_TRK2', 'class_TRK4', 'class_UND', 'num_roads', 'num_nodes',
                   'num_1degree', 'num_2degree', 'num_3degree', 'num_4degree',
                   'num_greater5degree']
node_df[road_attributes] += 1.0

# edit edge_df
print(np.sum(edge_df.isna()))
print(np.sum(edge_df==0))
print(edge_df.describe())
# drop self loops (110) and zero population areas (440)
non_zero_masks = np.logical_and(edge_df.od_duration > 0.000001, edge_df.total_pop_o > 0.000001) # od_duration and total_pop_o
non_zero_masks = np.logical_and(non_zero_masks, edge_df.total_pop_d > 0.000001)
non_zero_masks = np.logical_and(non_zero_masks, edge_df.avg_age_o > 0.000001)
non_zero_masks = np.logical_and(non_zero_masks, edge_df.avg_age_d > 0.000001)
non_zero_masks = np.logical_and(non_zero_masks, edge_df.bachelor_degree_percent_o > 0.000001)
non_zero_masks = np.logical_and(non_zero_masks, edge_df.bachelor_degree_percent_d > 0.000001)
non_zero_masks = np.logical_and(non_zero_masks, edge_df.poi_entropy_agg_o > 0.000001)
non_zero_masks = np.logical_and(non_zero_masks, edge_df.poi_entropy_agg_d > 0.000001)
edge_df = edge_df.loc[non_zero_masks, :]

# attach POI density (O and D) to edge_df
edge_df = edge_df.merge(node_df[['SA2_MAIN16', 'pop_density', 'poi_count_density', 'poi_count_agg_density', 'poi_entropy_density', 'poi_entropy_agg_density', 'job_density_1', 'job_density_2']],
                        left_on = 'O', right_on = 'SA2_MAIN16')
edge_df = edge_df.merge(node_df[['SA2_MAIN16', 'pop_density', 'poi_count_density', 'poi_count_agg_density', 'poi_entropy_density', 'poi_entropy_agg_density', 'job_density_1', 'job_density_2']],
                        left_on = 'D', right_on = 'SA2_MAIN16', suffixes = ['_o', '_d'])

# split edge_df into edge_consumption_df and edge_flow_df
# it is because the consumption data have a lot of missing info.
edge_flow_df = edge_df.loc[~edge_df.flow_agents.isna(), ]
edge_consumption_df = edge_df.loc[~edge_df.consumption_count_age_source.isna(), ]

#
edge_flow_df = edge_flow_df.drop(columns = ['consumption_count_age_source', 'consumption_amount_age_source',
                                            'consumption_count_mcc_source', 'consumption_amount_mcc_source'])
edge_flow_df.dropna(how='any', inplace=True)
edge_consumption_df.dropna(how='any', inplace=True)

# lift road attributes by one unit.
edge_flow_df[road_attributes] += 1.0
edge_consumption_df[road_attributes] += 1.0

# print dimensions
print("Final edge mobility flow dataframe dim is: ", edge_flow_df.shape)
print("Final edge consumption flow dataframe dim is: ", edge_consumption_df.shape)
print("Final node dataframe dim is: ", node_df.shape)

# check the final
print(np.sum(node_df.isna()))
print(np.sum(node_df == 0))
print(np.sum(edge_flow_df.isna()))
print(np.sum(edge_flow_df == 0))
print(np.sum(edge_consumption_df.isna()))
print(np.sum(edge_consumption_df == 0))

# save dataframes
node_df.to_pickle(processing_data_path+'node_df.pickle')
edge_flow_df.to_pickle(processing_data_path+'edge_flow_df.pickle')
edge_consumption_df.to_pickle(processing_data_path+'edge_consumption_df.pickle')
#endregion


text:
import numpy as np
import tensorflow as tf
from scipy.stats import beta
from matplotlib import pyplot as plt
from matplotlib.ticker import FormatStrFormatter

class StickBreakingProcess(object):

    def __init__(self):

        # set number of test samples to be used after fitting the distribution
        self.num_test_samples = int(1e5)

    def plot_marginal_dirichlet_fit(self, pi, alpha_target, d_kl, dist_name):
        """
        :param pi: samples of simplex random variables
        :param alpha_target: parameters of target Dirichlet
        :param d_kl: KL-Divergence objective (for plotting purposes)
        :param dist_name: name of distribution used to generate pi (for plotting purposes)
        :return: None
        """
        # get number of classes
        K = len(alpha_target)

        # get number of plots and determine plot arrangement
        num_plots = 1 + K
        num_rows = 3
        num_cols = int(np.ceil(num_plots / num_rows))

        # declare plot
        fig, ax = plt.subplots(num_rows, num_cols, figsize=(16, 10))
        ax = ax.reshape(-1)

        # plot learning curve
        ax[0].set_title('Training Objective: $D_{kl}(q || p)$')
        ax[0].plot(np.arange(1, 1 + len(d_kl)), d_kl, 'k-', lw=2, alpha=0.6, label='Dkl')
        ax[0].set_xlabel('Iteration')
        ax[0].set_ylabel('$D_{kl}(q || p)$')

        # take some samples from a true Dirichlet
        pi_true = np.random.dirichlet(alpha_target, self.num_test_samples)

        # compute useful term
        alpha_target_0 = np.sum(alpha_target)

        # loop over the Beta marginals
        for k in range(K):

            # subplot title
            ax[1 + k].set_title('Beta Marginal for $\\pi_{' + str(k + 1) + '}$')
            ax[1 + k].set_xlabel('$\\pi_{' + str(k + 1) + '}$')
            ax[1 + k].set_ylabel('Density')

            # set true Beta marginal evaluation points
            x = np.linspace(beta.ppf(0.01, alpha_target[k], alpha_target_0 - alpha_target[k]),
                            beta.ppf(0.99, alpha_target[k], alpha_target_0 - alpha_target[k]),
                            100)

            # plot true Beta marginal pdf
            ax[1 + k].plot(x, beta.pdf(x, alpha_target[k], alpha_target_0 - alpha_target[k]), 'k-', lw=2, label='pdf')

            # plot histogram densities of generated samples (proposed and true Dirichlet)
            hist = ax[1 + k].hist([pi[:, k], pi_true[:, k]],
                                  density=True,
                                  histtype='stepfilled',
                                  alpha=0.2,
                                  bins=100,
                                  label=(dist_name, 'Dirichlet'))

            # enable the legend
            ax[1 + k].legend()

            # limit y-axis appropriately for non-saturating Beta distributions
            ax[1 + k].set_ylim(0, 1.25 * max([np.max(x) for x in hist[0]]))

        # make it tight
        plt.tight_layout()

class KumaraswamyStickBreakingProcess(StickBreakingProcess):

    def __init__(self, dkl_taylor_order=10):
        """
        :param K: number of classes
        :param dkl_taylor_order: Taylor expansion approximation order for the KL-divergence computation
        """
        # initialize base class
        StickBreakingProcess.__init__(self)

        # save Taylor approximation order
        assert isinstance(dkl_taylor_order, int) and dkl_taylor_order >= 1
        self.M = dkl_taylor_order

    @staticmethod
    def __parameter_rank_check(alpha):
        """
        This function ensures any parameter vector has rank 2 (i.e. at least [1 x K])
        :param alpha: Dirichlet parameters
        :return: alpha: Dirichlet parameters
        """
        # if alpha is not already [batch size x K], ensure alpha is then shape [1 x K]
        if len(alpha.get_shape().as_list()) == 1:
            alpha = tf.expand_dims(alpha, axis=0)

        return alpha

    def __stick_break_parameters(self, alpha):
        """
        This function converts a Dirichlet alpha parameter vector to the Beta marginal parameters required by the
        stick-breaking Kumaraswamy-Dirichlet sampling procedure.
        :param alpha: Dirichlet parameters
        :return: a, b: marginal Beta distribution parameters
        """
        # enforce rank
        alpha = self.__parameter_rank_check(alpha)

        # get number of clusters
        K = alpha.shape.as_list()[-1]

        # compute stick-break parameters
        a = []
        b = []
        for k in range(K - 1):
            a.append(alpha[:, k])
            b.append(tf.reduce_sum(alpha[:, k+1:], axis=1))
        a = tf.stack(a, axis=1)
        b = tf.stack(b, axis=1)

        return a, b

    def sample(self, alpha, use_rand_perm=True):
        """
        This method implements the Kumarswamy-Dirichlet sample approximation via a stick-breaking process.
        :param alpha: Dirichlet parameters [batch size x K]
        :param use_rand_perm: whether to perform randomized permutation
        :return:
            pi: a Dirichlet sample approximation generated by string cutting Kumaraswamy samples
            i_perm: the sampled permutation indices
        """
        # if no input, return None
        if alpha is None:
            return [None] * 2

        # get number of clusters
        K = alpha.shape.as_list()[-1]

        # using randomized permutation
        if use_rand_perm:

            # sample permutation indices
            i_perm = tf.contrib.framework.argsort(tf.random_uniform(tf.shape(alpha), minval=0, maxval=1), axis=-1)
            i_perm = tf.reshape(i_perm, [-1, K])  # arg sort returns [None, None] despite documentation

            # apply permutation since the pi[K] has bias for sparse sampling
            alpha = tf.batch_gather(alpha, i_perm)

        # no random permutation
        else:

            i_perm = None

        # get stick-break parameters
        a, b = self.__stick_break_parameters(alpha)

        # sample uniform noise (setting minimum value > 0 resolves numerical stability issues)
        u = tf.random_uniform(tf.shape(a), minval=1e-4, maxval=1)

        # sample Kumaraswamy (Beta approximations) via inverse CDF sampling
        x = (1 - (1 - u) ** (1 / b)) ** (1 / a)

        # convert to a Dirichlet sample approximation via Beta string cutting method
        pi = [x[:, 0]]
        for j in range(1, K - 1):
            pi.append((1 - tf.reduce_sum(tf.stack(pi, axis=-1), axis=-1)) * x[:, j])
        pi.append(1 - tf.reduce_sum(tf.stack(pi, axis=-1), axis=-1))
        pi = tf.stack(pi, axis=1)

        # using randomized permutation
        if use_rand_perm:

            # construct inverse permutation indices
            i_inv_perm = tf.contrib.framework.argsort(i_perm, axis=-1)
            i_inv_perm = tf.reshape(i_inv_perm, [-1, K])  # once again, arg sort returns [None, None]

            # invert permutation
            pi = tf.batch_gather(pi, i_inv_perm)

        return pi, i_perm

    def kl_divergence(self, alpha, alpha_prior, i_perm=None, wrt='Dirichlet-Marginals'):
        """
        Computes the KL divergence between the Kumaraswamy q distributions and the Dirichlet prior's Beta marginals.
        :param alpha: posterior approximation Dirichlet parameters
        :param alpha_prior: prior Dirichlet parameters
        :param i_perm: random permutation indices used during sampling procedure
        :param wrt: that which the KL divergence is with respect to, either Dirichlet marginal or Beta stick breaks
        :return: KL divergence of marginal Beta distributions of shape [batch size x K]
        """
        assert wrt in {'Dirichlet-Marginals', 'Beta-Sticks'}

        # apply permutation if one was provided
        if i_perm is not None:
            alpha_prior = self.__parameter_rank_check(alpha_prior)
            alpha_prior = tf.tile(alpha_prior, tf.stack((tf.shape(alpha)[0], 1)))
            alpha = tf.batch_gather(alpha, i_perm)
            alpha_prior = tf.batch_gather(alpha_prior, i_perm)

        # take KL divergence w.r.t. to the Dirichlet's marginal Betas
        if wrt == 'Dirichlet-Marginals':

            # compute marginal q(pi; a', b') approximation parameters
            a_prime = self.__parameter_rank_check(alpha)
            b_prime = tf.reduce_sum(a_prime, axis=1, keepdims=True) - a_prime

            # compute marginal p(pi; a, b) prior parameters
            a_prior = self.__parameter_rank_check(alpha_prior)
            b_prior = tf.reduce_sum(a_prior, axis=1, keepdims=True) - a_prior

        # take KL divergence w.r.t. to the stick-breaking marginal Betas
        else:

            # compute marginal q(pi; a', b') approximation parameters
            a_prime, b_prime = self.__stick_break_parameters(alpha)

            # compute marginal p(pi; a, b) prior parameters
            a_prior, b_prior = self.__stick_break_parameters(alpha_prior)

        # KL-Divergence
        kl = (a_prime - a_prior) / a_prime * (-np.euler_gamma - tf.digamma(b_prime) - 1 / b_prime) \
            + (tf.log(a_prime * b_prime)) \
            + (tf.lbeta(tf.stack((a_prior, b_prior), axis=-1))) \
            - (b_prime - 1) / b_prime
        for m in range(1, self.M + 1):
            B = tf.exp(tf.lbeta(tf.concat((tf.expand_dims(m / a_prime, axis=-1), tf.expand_dims(b_prime, axis=-1)), axis=-1)))
            kl += (b_prior - 1) * b_prime / (m + a_prime * b_prime) * B

        # sum over the dimensions
        kl = tf.reduce_sum(kl, axis=1)

        return kl

    def unit_tests(self, K=10):
        """
        Runs some unit tests on both the sampling procedure and the KL divergence computation.
        :param K: number of classes to run tests over
        :return: None
        """
        # begin test session
        tf.reset_default_graph()
        cfg = tf.ConfigProto()
        cfg.gpu_options.allow_growth = True
        with tf.Session(config=cfg) as sess:

            # declare parameter placeholder
            alpha_ph = tf.placeholder(dtype=tf.float32, shape=[self.num_test_samples, K])

            # sample some random values
            alpha = 10 * np.random.ranf([self.num_test_samples, K])

            # ensure probabilities sum to 1
            pi = sess.run(self.sample(alpha_ph)[0], feed_dict={alpha_ph: alpha})
            assert (np.abs(np.sum(pi, axis=1) - 1) <= 1e-6).all()

            # set an appropriate prior
            alpha_prior = tf.constant(np.ones(K) / K, dtype=tf.float32)

            # ensure KL divergences are non-negative
            kl = sess.run(self.kl_divergence(alpha_ph, alpha_prior), feed_dict={alpha_ph: alpha})
            assert (kl >= 0).all()

            # test batch gather functionality
            x = np.array([[0, 1, 2], [0, 2, 1],
                          [1, 0, 2], [1, 2, 0],
                          [2, 0, 1], [2, 1, 0]])
            x_sort = np.zeros(x.shape)
            for i in range(x.shape[0]):
                x_sort[i] = x[i, x[i]]
            x_test = sess.run(tf.batch_gather(tf.constant(x, dtype=tf.int32), tf.constant(x, dtype=tf.int32)))
            assert (x_sort == x_test).all()

        # print success
        print('Kumaraswamy-Dirichlet Sampler unit tests passed!')

    def fit(self, alpha_target, use_rand_perm, kl_wrt, learning_rate=5e-3, num_epochs=int(5e4)):
        """
        This is a test method that minimizes the KL divergence between the Kumarswamy approximations to the target
        Dirichlet's Beta marginals. It will plot the true Beta marginal PDF over histograms for exact marginalized Beta
        samples and samples generated by the Kumarswamy-Dirichlet sampler.
        :param alpha_target: target Dirichlet parameters
        :param learning_rate: learning rate applied to KL-divergence gradient
        :param num_epochs: number of training epochs
        :param use_rand_perm: whether to use random permutation matrices during sampling and KL divergence computation
        :param kl_wrt: that which the KL divergence is with respect to, either Dirichlet marginal or Beta stick breaks
        :return: samples in the simplex
        """
        # get number of classes
        K = len(alpha_target)

        # begin test session
        tf.reset_default_graph()
        cfg = tf.ConfigProto()
        cfg.gpu_options.allow_growth = True
        with tf.Session(config=cfg) as sess:

            # declare parameters we wish to fit
            alpha = tf.get_variable('alpha',
                                    shape=[1, K],
                                    dtype=tf.float32,
                                    initializer=tf.random_uniform_initializer(minval=0, maxval=15))

            # define loss operation
            if use_rand_perm:
                kl = self.kl_divergence(alpha, tf.constant(alpha_target, dtype=tf.float32),
                                        i_perm=self.sample(alpha)[1],
                                        wrt=kl_wrt)
            else:
                kl = self.kl_divergence(alpha, tf.constant(alpha_target, dtype=tf.float32),
                                        i_perm=None,
                                        wrt=kl_wrt)
            loss_op = tf.reduce_sum(kl)

            # configure optimizer
            optimizer = tf.train.GradientDescentOptimizer(learning_rate)

            # configure training operation
            train_op = tf.contrib.layers.optimize_loss(loss=loss_op,
                                                       global_step=tf.Variable(0, name='global_step', trainable=False),
                                                       learning_rate=learning_rate,
                                                       optimizer=optimizer)

            # run initialization
            sess.run(tf.global_variables_initializer())

            # train it
            loss = np.zeros(num_epochs)
            for t in range(num_epochs):
                _, loss[t] = sess.run([train_op, loss_op])
            print('Final DKL = {:.2f}'.format(loss[-1]))
            print('Alpha = ', sess.run(alpha))

            # take some samples
            pi = sess.run(self.sample(tf.tile(alpha, tf.constant([self.num_test_samples, 1])),
                                      use_rand_perm=use_rand_perm)[0])

            # plot results
            self.plot_marginal_dirichlet_fit(pi, alpha_target, loss, 'KDS')

        return pi

class BetaStickBreakingProcess(StickBreakingProcess):

    def __init__(self):
        # initialize base class
        StickBreakingProcess.__init__(self)

    @staticmethod
    def sample(alpha, n):
        """
        This method implements the Dirichlet samples via a stick-breaking process using Beta distributions.
        :param alpha: Dirichlet parameters
        :return: pi: a Dirichlet sample
        """
        # squeeze any extra alpha dimensions
        alpha = np.squeeze(alpha)

        # initialize x
        x = np.zeros([n, len(alpha)])
        pi = np.zeros([n, len(alpha)])

        # take initial sample
        x[:, 0] = np.random.beta(alpha[0], np.sum(alpha[1:]), n)

        # convert to a Dirichlet sample approximation via string cutting
        pi[:, 0] = x[:, 0]
        for j in range(1, len(alpha) - 1):
            x[:, j] = np.random.beta(alpha[j], np.sum(alpha[j + 1:]), n)
            pi[:, j] = ((1 - np.sum(pi[:, :j], axis=-1)) * x[:, j])
        pi[:, -1] = (1 - np.sum(pi[:, :-1], axis=-1))

        # make sure they sum to 1
        assert (np.sum(pi, axis=1) == 1).all()

        return pi

    @staticmethod
    def kl_divergence(alpha, alpha_prior):
        """
        Computes the KL divergence between two Dirichlet distributions
        :param alpha: posterior approximation Dirichlet parameters
        :param alpha_prior: prior Dirichlet parameters
        :return: KL divergence
        """
        # compute convenient terms
        alpha_0 = tf.reduce_sum(alpha, axis=-1, keepdims=True)
        beta_0 = tf.reduce_sum(alpha_prior, axis=-1, keepdims=True)

        # compute KL(q || p)
        kl = \
            tf.lgamma(alpha_0) - \
            tf.reduce_sum(tf.lgamma(alpha), axis=-1) - \
            tf.lgamma(beta_0) + \
            tf.reduce_sum(tf.lgamma(alpha_prior)) + \
            tf.reduce_sum((alpha - alpha_prior) * (tf.digamma(alpha) - tf.digamma(alpha_0)), axis=-1)

        return kl

    def fit(self, alpha_target, learning_rate=1e-2, num_epochs=int(1e4)):
        """
        :param alpha_target: target Dirichlet parameters
        :param learning_rate: learning rate applied to KL-divergence gradient
        :param num_epochs: number of training epochs
        :return: None
        """
        # get number of classes
        K = len(alpha_target)

        # begin test session
        tf.reset_default_graph()
        cfg = tf.ConfigProto()
        cfg.gpu_options.allow_growth = True
        with tf.Session(config=cfg) as sess:

            # declare parameters we wish to fit
            alpha = tf.get_variable('alpha',
                                    shape=[1, K],
                                    dtype=tf.float32,
                                    initializer=tf.random_uniform_initializer(minval=0, maxval=15))

            # define loss operation
            loss_op = tf.reduce_sum(self.kl_divergence(alpha, tf.constant(alpha_target, dtype=tf.float32)))

            # configure optimizer
            optimizer = tf.train.AdamOptimizer(learning_rate)

            # configure training operation
            train_op = tf.contrib.layers.optimize_loss(loss=loss_op,
                                                       global_step=tf.Variable(0, name='global_step', trainable=False),
                                                       learning_rate=learning_rate,
                                                       optimizer=optimizer)

            # run initialization
            sess.run(tf.global_variables_initializer())

            # train it
            loss = np.zeros(num_epochs)
            for t in range(num_epochs):
                _, loss[t] = sess.run([train_op, loss_op])
            print('Final DKL = {:.2f}'.format(loss[-1]))

            # take some samples
            alpha_fit = sess.run(alpha)
            print(alpha_fit)
            pi = self.sample(alpha_fit, self.num_test_samples)

            # plot results
            self.plot_marginal_dirichlet_fit(pi, alpha_target, loss, 'KBS')

def mv_kumaraswamy_ordering_impact(pi_no_perm, pi_random_perm, alpha_target):
    """
    Plots MV Kumaraswamy's ordering's impact on bias (i.e. plots the marginals and shows bias on last dimension)
    :param pi_no_perm: samples that did not use random permutations
    :param pi_random_perm: samples that used random permutations
    :param alpha_target: Dirichlet target parameters
    :return: None
    """
    # get number of classes
    K = len(alpha_target)

    # declare plot
    fig, ax = plt.subplots(1, K, figsize=(8, 2))
    ax = ax.reshape(-1)

    # take some samples from a true Dirichlet
    pi_true = np.random.dirichlet(alpha_target, int(1e6))

    # compute useful term
    alpha_target_0 = np.sum(alpha_target)

    # loop over the Beta marginals
    for k in range(K):

        # enable y label for first axis only
        if K == 0:
            ax[k].set_ylabel('Density')

        # set true Beta marginal evaluation points
        x = np.linspace(beta.ppf(0.01, alpha_target[k], alpha_target_0 - alpha_target[k]),
                        beta.ppf(0.99, alpha_target[k], alpha_target_0 - alpha_target[k]),
                        int(1e5))

        # plot true Beta marginal pdf
        ax[k].plot(x, beta.pdf(x, alpha_target[k], alpha_target_0 - alpha_target[k]), 'k:', lw=2, label='PDF')

        # plot histogram densities of generated samples (proposed and true Dirichlet)
        hist = ax[k].hist([pi_true[:, k], pi_no_perm[:, k], pi_random_perm[:, k]],
                          density=True,
                          histtype='bar',
                          alpha=0.2,
                          bins=100,
                          label=('Dirichlet',
                                 '$f_{12345}$',
                                 '$E[f]$'))

        # limit axes appropriately
        ax[k].set_xlim(0, 0.05)
        ax[k].set_ylim(0, 55)

        # configure ticks
        ax[k].set_xticks([0, 0.03])
        if k > 0:
            ax[k].set_yticks([])
        ax[k].xaxis.set_major_formatter(FormatStrFormatter('%.2f'))
        ax[k].yaxis.set_major_formatter(FormatStrFormatter('%.0f'))

        # enable the legend
        if k == int(K / 2):
            ax[k].legend(ncol=4, bbox_to_anchor=(2, -0.125))

    # make it tight
    plt.subplots_adjust(left=0.04, bottom=0.3, right=0.98, top=0.98, wspace=0.05, hspace=0)

if __name__ == '__main__':

    # declare Kumaraswamy stick-breaking sampler and run unit tests
    ksb = KumaraswamyStickBreakingProcess()
    ksb.unit_tests()

    # demonstrate effect ordering has on approximating sparsity-inducing Dirichlet
    alpha = np.ones(5) / 5
    pi_no_perm = ksb.fit(alpha, use_rand_perm=False, kl_wrt='Dirichlet-Marginals')
    pi_random_perm = ksb.fit(alpha, use_rand_perm=True, kl_wrt='Dirichlet-Marginals')
    mv_kumaraswamy_ordering_impact(pi_no_perm, pi_random_perm, alpha)
    plt.show()

    # declare Beta stick-breaking sampler (samples a Dirichlet random variable)
    bsb = BetaStickBreakingProcess()

    # declare some Dirichlet parameters to fit
    alphas = [np.ones(5),
              np.ones(5) / 5,
              1 + np.arange(5),
              np.flip(1 + np.arange(5)),
              (1 + np.arange(5)) * 2,
              np.flip((1 + np.arange(5)) * 2)]

    # loop over the parameters
    for a in alphas:

        # fit Kumaraswamy-Dirichlet sample under various methods
        ksb.fit(a, use_rand_perm=False, kl_wrt='Dirichlet-Marginals')  # not accurate w.r.t the model
        ksb.fit(a, use_rand_perm=True, kl_wrt='Dirichlet-Marginals')  # same as above--perm. does not affect KL wrt Dir.
        ksb.fit(a, use_rand_perm=False, kl_wrt='Beta-Sticks')  # without random permutation, we incur bias
        ksb.fit(a, use_rand_perm=True, kl_wrt='Beta-Sticks')   # truest to model--but stochastic (breaks DKL convexity)

        # fit Beta-Dirichlet sampler
        bsb.fit(a)

        # show results
        plt.show()

max_stars_repo_name: astirn/KumaraswamyVAE

max_issues_repo_name: astirn/KumaraswamyVAE

max_forks_repo_name: astirn/KumaraswamyVAE

reason: import numpy,from scipy


text:
from angela2 import *
from angela2.utils.logger import catch_logs
from numpy.testing import assert_equal, assert_array_equal
import pytest

@pytest.mark.codegen_independent
def test_clock_attributes():
    clock = Clock(dt=1*ms)
    assert_array_equal(clock.t, 0*second)
    assert_array_equal(clock.timestep, 0)
    assert_array_equal(clock.dt, 1*ms)

@pytest.mark.codegen_independent
def test_clock_dt_change():
    clock = Clock(dt=1*ms)
    # at time 0s, all dt changes should be allowed
    clock.dt = 0.75*ms
    clock._set_t_update_dt()
    clock.dt = 2.5*ms
    clock._set_t_update_dt()
    clock.dt = 1*ms
    clock._set_t_update_dt()

    # at 0.1ms only changes that are still representable as an integer of the
    # current time 1s are allowed
    clock.dt = 0.1*ms
    clock._set_t_update_dt()

    clock.dt = 0.05*ms
    clock._set_t_update_dt(target_t=0.1*ms)
    clock.dt = 0.1*ms
    clock._set_t_update_dt(target_t=0.1*ms)
    clock.dt = 0.3*ms
    with pytest.raises(ValueError):
        clock._set_t_update_dt(target_t=0.1*ms)

@pytest.mark.codegen_independent
def test_defaultclock():
    defaultclock.dt = 1*ms
    assert_equal(defaultclock.dt, 1*ms)
    assert defaultclock.name == 'defaultclock'

@pytest.mark.codegen_independent
def test_set_interval_warning():
    clock = Clock(dt=0.1*ms)
    with catch_logs() as logs:
        clock.set_interval(0*second, 1000*second)  # no problem
    assert len(logs) == 0
    with catch_logs() as logs:
        clock.set_interval(0*second, 10000000000*second)  # too long
    assert len(logs) == 1
    assert logs[0][1].endswith('many_timesteps')

if __name__ == '__main__':
    test_clock_attributes()
    restore_initial_state()
    test_clock_dt_change()
    restore_initial_state()
    test_defaultclock()
    test_set_interval_warning()

reason: from numpy


text:
'''
Ensemble multiple model:
1. concat
--------------------------------------
Author: Yang Qian
Email: yqian@aibee.com
'''
import argparse
import logging
import sys
import os
import atexit
import bisect

import numpy as np
import cv2
import torch
import tqdm
from torch.backends import cudnn
from collections import deque
from collections import OrderedDict

sys.path.append('.')

from fastreid.evaluation import evaluate_rank, ReidEvaluator, print_csv_format
from fastreid.config import get_cfg
from fastreid.utils.logger import setup_logger
from fastreid.data import build_reid_test_loader
from fastreid.utils.visualizer import Visualizer
from fastreid.utils.events import TensorboardXWriter
from fastreid.modeling.meta_arch import build_model
from fastreid.utils.checkpoint import Checkpointer
from fastreid.engine import DefaultPredictor

cudnn.benchmark = True
logger = logging.getLogger('fastreid.ensembel')

def setup_cfg(config_file, opts):
    # load config from file and command-line arguments
    cfg = get_cfg()
    cfg.merge_from_file(config_file)
    cfg.merge_from_list(opts)
    cfg.freeze()
    return cfg

def get_parser():
    parser = argparse.ArgumentParser(description="fastreid model ensembel")
    parser.add_argument("--config-file", default=["configs/Aibee/Store-CTF-body.yml"], metavar="FILE", nargs='+', help="path to config file")
    parser.add_argument(
        "--resume",
        action="store_true",
        help="whether to attempt to resume from the checkpoint directory",
    )

    parser.add_argument(
        "opts",
        help="Modify config options using the command-line",
        default=None,
        nargs=argparse.REMAINDER,
    )
    return parser

class FeatureExtraction(object):
    def __init__(self, cfgs):
        """
        Args:
            cfg (CfgNode):
        """
        self.cfgs = cfgs
        self.predictors = []

        for cfg in cfgs:
            self.predictors.append(DefaultPredictor(cfg))

    def run_on_loader(self, data_loader):
        for batch in data_loader:
            predictions = []
            for predictor in self.predictors:
                predictions.append(predictor(batch["images"]))
            yield torch.cat(predictions, dim=-1), batch

if __name__ == "__main__":
    args = get_parser().parse_args()
    logger = setup_logger()
    cfgs = []
    for config_file in args.config_file:
        cfg = setup_cfg(config_file, args.opts)
        cfgs.append(cfg)
    results = OrderedDict()
    for dataset_name in cfgs[0].DATASETS.TESTS:
        test_loader, num_query = build_reid_test_loader(cfgs[0], dataset_name)
        evaluator = ReidEvaluator(cfgs[0], num_query)
        feat_extract = FeatureExtraction(cfgs)
        for (feat, batch) in tqdm.tqdm(feat_extract.run_on_loader(test_loader), total=len(test_loader)):
            evaluator.process(batch, feat)
        result = evaluator.evaluate()
        results[dataset_name] = result
    print_csv_format(results)


text:
#!/usr/bin/env python3

# apply_surprise_metric.py

import sys, os, csv, random
import numpy as np
from collections import Counter

date = sys.argv[1]
novelname = sys.argv[2]

# USAGE:
# python3 apply_normalized_surprise_metric.py 1840 leguinlefthand.txt

# The program looks for a surprise-metric file in the surprisemetrics subfolder,
# using "date" (above) to generate the name of the file.

# It also looks for "novelname" in a sibling folder called "rawtexts."

# Then it applies the surprise-metric to the novel, after first breaking
# the novel into ~250-word chunks. It uppercases the 2000 words that have the
# biggest positive effect on surprise (i.e., words that are more likely
# in the later of the two periods that were contrasted to create the
# the surprise metric).

def onlyalpha(word):
    only = ''
    for char in word:
        if char.isalpha():
            only = only + char
    return only

inpath = '../rawtexts/' + novelname

with open(inpath, encoding = 'utf-8') as f:
    lines = f.readlines()

coefs = dict()
with open('surprisemetrics/surprise_in_' + str(date) + '.csv', encoding = 'utf-8') as f:
    reader = csv.DictReader(f)
    for row in reader:
        coefs[row['word']] = float(row['coef'])

pages = []
ctr = 0
for l in lines:
    if l.startswith('“Stranger In'):
        continue

    if ctr == 0:
        page = []
    page.append(l)
    wordct = len(l.split())
    ctr += wordct
    if ctr > 250:
        pages.append(page)
        ctr = 0

# first let's get total wordcounts
wordcounts = Counter()

for p in pages:

    for l in p:
        l = l.replace('-', ' ')
        l = l.replace('—', ' ')
        l = l.replace(',', ' ')
        words = l.split()
        for origw in words:
            w = onlyalpha(origw.lower())
            if w in coefs:
                wordcounts[w] += 1

significant = []

for word, coef in coefs.items():
    if coef > 0:
        significant.append((coef, word))

significant.sort(reverse = True)

significant = set([x[1] for x in significant[0: 2000]])

pageprobs = []
transpages = []

for p in pages:
    wordcounts = Counter()
    total = 0

    translated = []

    for l in p:
        l = l.replace('-', ' ')
        l = l.replace('—', ' ')
        l = l.replace(',', ' ')
        words = l.split()
        for origw in words:
            w = onlyalpha(origw.lower())
            if w in coefs:
                wordcounts[w] += 1
                if w in significant:
                    translated.append(origw.upper())
                else:
                    translated.append(origw)
            else:
                translated.append(origw)

            total += 1

    prob = 0

    for word, count in wordcounts.items():
        freq = count / total
        prob += (freq * coefs[word])

    pageprobs.append(prob)
    transpages.append(' '.join(translated))

tuplepages = [x for x in zip(pageprobs, transpages)]

tuplepages.sort(reverse = True)

with open('surprise_' + novelname, mode = 'w', encoding = 'utf-8') as f:
    for p, t in tuplepages:
        f.write('SURPRISE: ' + str(p) + '\n')
        f.write('\n')
        f.write(t)
        f.write('\n\n')

max_stars_repo_name: tedunderwood/measureperspective

max_issues_repo_name: tedunderwood/measureperspective

max_forks_repo_name: tedunderwood/measureperspective


text:
from importlib.resources import open_binary
import numpy as np
import struct

def get_mapping():
    '''
    Generates a mapping of Myanmar Unicode characters to their corresponding
    indices in the parameter array.

    Returns
    -------
    dict
        A mapping from Myanmar Unicode characters to indices.
    '''
    def get_chars(start, end):
        return [chr(char) for char in range(ord(start), ord(end) + 1)]

    chars = (get_chars('\u1000', '\u103F')
            + get_chars('\u104A', '\u109F')
            + get_chars('\uAA60', '\uAA7F')
            + get_chars('\uA9E0', '\uA9FF')
            + get_chars('\u2000', '\u200B'))

    return {char: i + 1 for i, char in enumerate(chars)}

def load_params():
    '''
    Loads parameters as a 2d array, which are log likelihood ratios of
    Unicode to Zawgyi.

    Returns
    -------
    numpy.ndarray
        Parameters as a 2d array.
    '''
    def read_char_array(f, size):
        return struct.unpack(f'{size}s', f.read(size))[0].decode('utf-8')

    def read_float(f):
        return struct.unpack('>f', f.read(4))[0]

    def read_int(f):
        return struct.unpack('>i', f.read(4))[0]

    def read_short(f):
        return struct.unpack('>h', f.read(2))[0]

    with open_binary('myanmartools.resources', 'zawgyiUnicodeModel.dat') as f:
        # check signature
        uzmodel_tag = read_char_array(f, 8)
        if uzmodel_tag != 'UZMODEL ':
            raise IOError('incorrect uzmodel_tag')
        uzmodel_version = read_int(f)
        if uzmodel_version == 1:
            ssv = 0
        elif uzmodel_version == 2:
            ssv = read_int(f)
        else:
            raise IOError('incorrect uzmodel_version')
        bmarkov_tag = read_char_array(f, 8)
        if bmarkov_tag != 'BMARKOV ':
            raise IOError('incorrect bmarkov_tag')
        bmarkov_version = read_int(f)
        if bmarkov_version != 0:
            raise IOError('incorrect bmarkov_version')

        # read params
        size = read_short(f)
        params = np.empty((size, size))
        for row in range(size):
            count = read_short(f)
            if count != 0:
                params[row] = read_float(f)
                for i in range(count):
                    col = read_short(f)
                    params[row, col] = read_float(f)
            else:
                params[row] = 0

        return params

max_stars_repo_name: sven-oly/myanmar-tools

max_issues_repo_name: sven-oly/myanmar-tools

max_forks_repo_name: sven-oly/myanmar-tools


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri May  8 15:42:06 2020

@author: degerli
"""
import os
import argparse
import numpy as np
import imageio

def convertNumpy(path):
    data_list = sorted(os.listdir(path))
    data = []
    for i in data_list:
        #save data
        img = imageio.imread(path + i)
        img_resized_float = img.astype('float32')
        data.append(np.expand_dims(np.array(img_resized_float), axis = 2))

    data = np.array(data)
    return data

ap = argparse.ArgumentParser()
ap.add_argument('--dir_early', required=True, help="path to early COVID-19 images")
ap.add_argument('--dir_mild', required=True, help="path to mild COVID-19 images")
ap.add_argument('--dir_control', required=True, help="path to control group images")
args = vars(ap.parse_args())

np.save(('data/covid_early.npy'), convertNumpy(args['dir_early']))
np.save(('data/covid_mild.npy'), convertNumpy(args['dir_mild']))
np.save(('data/normal.npy'), convertNumpy(args['dir_control']))


text:
#coding=utf-8
"""Module for processing data

The function of this Module is served for processing data.

"""
import pandas as pd
import numpy as np

def cut_groups(data, col, cutoffs):
    """Cut data into subsets according to cutoffs

    Parameters
    ----------
    data : pandas.DataFrame
        Data to split.
    col : str
        Name of column in data to compare with.
    cutoffs : list(int)
        List of cutoffs, like as [min-value, 30, 60, max-value].

    Returns
    -------
    list(pandas.DataFrame)
        List of sub-data as DataFrame.

    Examples
    --------
    >>> cut_groups(data, "X", [0, 0.4, 0.6, 1.0])
    [pandas.DataFrame, pandas.DataFrame, pandas.DataFrame]
    """
    res = []
    N = len(cutoffs)
    for i in range(N - 1):
        if i == N - 2:
            df = data[(data[col] >= cutoffs[i]) & (data[col] <= cutoffs[i+1])]
        else:
            df = data[(data[col] >= cutoffs[i]) & (data[col] < cutoffs[i+1])]
        res.append(df)
    return res

def prepare_surv(x, label):
    """Prepare data for survival analyze(Deep Surival).

    Parameters
    ----------
    x : numpy.array
        Two-dimension array indicating variables.
    label : dict
        Contain 'e', 't'.

        examples as {'e': np.array, 't': np.array}.

    Returns
    -------
    tuple
        Sorted (x, label) tuple of survival data.

    Examples
    --------
    >>> prepare_surv(data[x_cols].values, {'e': data['e'].values, 't': data['t'].values})
    """
    if isinstance(label, dict):
       e, t = label['e'], label['t']

    # Sort Training Data for Accurate Likelihood
    # sort array using pandas.DataFrame(According to DESC 't' and ASC 'e')
    df1 = pd.DataFrame({'t': t, 'e': e})
    df1.sort_values(['t', 'e'], ascending=[False, True], inplace=True)
    sort_idx = list(df1.index)
    x = x[sort_idx]
    e = e[sort_idx]
    t = t[sort_idx]

    return x, {'e': e, 't': t}

def parse_surv(x, label):
    """Parse raw-data for survival analyze(Deep Surival).

    Parameters
    ----------
    x : np.array
        two-dimension array indicating variables.
    label : dict
        Contain 'e', 't'.

        examples as {'e': np.array, 't': np.array}.

    Returns
    -------
    tuple
        Sorted (x, e, t) tuple, index of people who is failure or at risk, and type of ties.

    Examples
    --------
    >>> parse_surv(data[x_cols].values, {'e': data['e'].values, 't': data['t'].values})
    """
    # sort data by t
    x, label = prepare_surv(x, label)
    e, t = label['e'], label['t']

    failures = {}
    atrisk = {}
    n, cnt = 0, 0

    for i in range(len(e)):
        if e[i]:
            if t[i] not in failures:
                failures[t[i]] = [i]
                n += 1
            else:
                # ties occured
                cnt += 1
                failures[t[i]].append(i)

            if t[i] not in atrisk:
                atrisk[t[i]] = []
                for j in range(0, i+1):
                    atrisk[t[i]].append(j)
            else:
                atrisk[t[i]].append(i)
    # when ties occured frequently
    if cnt >= n / 2:
        ties = 'efron'
    elif cnt > 0:
        ties = 'breslow'
    else:
        ties = 'noties'

    return x, e, t, failures, atrisk, ties


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Author: Yuki Furuta <furushchev@jsk.imi.i.u-tokyo.ac.jp>

import os
import numpy as np
from .. import utils

url = "http://www.cs.toronto.edu/~nitish/unsupervised_video/mnist_test_seq.npy"

def get_moving_mnist():
    return utils.cache_or_load_file(
        "moving_mnist", url, os.path.basename(url), np.load)

if __name__ == '__main__':
    data = get_moving_mnist()
    print(data.shape)

max_stars_repo_name: otsubo/chainervr

max_issues_repo_name: otsubo/chainervr

max_forks_repo_name: otsubo/chainervr