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

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
import tensorflow as tf

def reflection(data, decision):
  """Conditionally reflects the data in XYZ.

  Args:
    data: input tensor, shape: [..], z, y, x, c
    decision: boolean tensor, shape 3, indicating on which spatial dimensions
       to apply the reflection (x, y, z)

  Returns:
    TF op to conditionally apply reflection.
  """
  with tf.compat.v1.name_scope('augment_reflection'):
    rank = data.get_shape().ndims
    spatial_dims = tf.constant([rank - 2, rank - 3, rank - 4])
    selected_dims = tf.boolean_mask(tensor=spatial_dims, mask=decision)
    return tf.reverse(data, selected_dims)

def xy_transpose(data, decision):
  """Conditionally transposes the X and Y axes of a tensor.

  Args:
    data: input tensor, shape: [..], y, x, c.
    decision: boolean scalar indicating whether to apply the transposition

  Returns:
    TF op to conditionally apply XY transposition.
  """
  with tf.compat.v1.name_scope('augment_xy_transpose'):
    rank = data.get_shape().ndims
    perm = range(rank)
    perm[rank - 3], perm[rank - 2] = perm[rank - 2], perm[rank - 3]
    return tf.cond(pred=decision,
                   true_fn=lambda: tf.transpose(a=data, perm=perm),
                   false_fn=lambda: data)

def permute_axes(x, permutation, permutable_axes):
  """Permutes the axes of `x` using the specified permutation.

  All axes not in `permutable_axes` must be identity mapped by `permutation`:
    `permutation[i] = i   if i not in permutable_axes`

  Args:
    x: Tensor of rank N to permute.
    permutation: Rank 1 tensor of shape `[N]`.
    permutable_axes: List (not Tensor) of axes that may be permuted.

  Returns:
    Permuted Tensor obtained by calling:
      `tf.transpose(x, permutation)`
    but with additional static shape information due to the restriction imposed
    by `permutable_axes`.
  """
  x = tf.convert_to_tensor(value=x)
  shape = x.shape.as_list()
  result = tf.transpose(a=x, perm=permutation)
  # Transpose loses shape information because it does not know that the
  # permutation is restricted to the permutable_axes, but we can restore
  # some or all of it.

  # If the static dimensions of all of the permutable axes are not the
  # same, then set them all to None.  Otherwise, they can stay as is.
  if len({shape[i] for i in permutable_axes}) > 1:
    for i in permutable_axes:
      shape[i] = None

  result.set_shape(shape)
  return result

class PermuteAndReflect(object):
  """Class for performing random permutation and reflection of axes.

  Consructing an instance of this class adds tensors to the default Tensorflow
  graph representing a randomly sampled permutation of `permutable_axes` and
  randomly sampled reflection decisions for `reflectable_axes`.

  Calling an instance of this class as a function (i.e. invoking the `__call__`
  method) applies the sampled transformations to a specified `Tensor`.

  Attributes:
    rank: The rank of Tensor that can be transformed.
    permutable_axes: 1-D int32 numpy array specifying the axes that may be
      permuted.
    reflectable_axes: 1-D int32 numpy array specifying the axes that may be
      reflected.
    reflect_decisions: bool Tensor of shape `[len(reflectable_axes)]` containing
      the sampled reflection decision for each axis in `reflectable_axes`.
    reflected_axes: Rank 1 int32 Tensor specifying the axes to be reflected
      (i.e. corresponding to a True value in `reflect_decisions`).
    permutation: int32 Tensor of shape `[len(permutable_axes)]` containing the
      sampled permutation of `permutable_axes`.
    full_permutation: int32 Tensor of shape `[rank]` that extends `permutation`
      to be a permutation of all `rank` axes, where axes not in
      `permutable_axes` are identity mapped.
  """

  def __init__(self, rank, permutable_axes, reflectable_axes,
               permutation_seed=None, reflection_seed=None):
    """Initializes the transformation nodes.

    Args:
      rank: The rank of the Tensor to be transformed.
      permutable_axes: The list (not a Tensor) of axes to be permuted.
      reflectable_axes: The list (not a Tensor) of axes to be reflected.
      permutation_seed: Optional integer.  Seed value to use for sampling axes
        permutation.
      reflection_seed: Optional integer.  Seed value to use for sampling
        reflection decisions.
    Raises:
      ValueError: if arguments are invalid.
    """
    self.rank = rank
    if len(set(permutable_axes)) != len(permutable_axes):
      raise ValueError('permutable_axes must not contain duplicates')
    if len(set(reflectable_axes)) != len(reflectable_axes):
      raise ValueError('reflectable_axes must not contain duplicates')
    if not all(0 <= x < rank for x in permutable_axes):
      raise ValueError('permutable_axes must be a subset of [0, rank-1].')
    if not all(0 <= x < rank for x in reflectable_axes):
      raise ValueError('reflectable_axes must be a subset of [0, rank-1].')
    # Cast to int32 to ensure the proper tensorflow types are used.
    self.permutable_axes = np.array(permutable_axes, dtype=np.int32)
    self.reflectable_axes = np.array(reflectable_axes, dtype=np.int32)

    if self.reflectable_axes.size > 0:
      self.reflect_decisions = tf.random.uniform([len(self.reflectable_axes)],
                                                 seed=reflection_seed) > 0.5
      self.reflected_axes = tf.boolean_mask(tensor=self.reflectable_axes,
                                            mask=self.reflect_decisions)

    if self.permutable_axes.size > 0:
      self.permutation = tf.random.shuffle(self.permutable_axes,
                                           seed=permutation_seed)
      # full_permutation must be a list rather than an np.array of int32 because
      # some elements are set to be tensors below.
      full_permutation = [np.int32(x) for x in range(rank)]
      for i, d in enumerate(self.permutable_axes):
        full_permutation[d] = self.permutation[i]
      self.full_permutation = tf.stack(full_permutation)

  def __call__(self, x):
    """Applies the sampled permutation and reflection to `x`.

    Args:
      x: A Tensor of rank `self.rank`.

    Returns:
      The transformed Tensor, retaining as much static shape information as
      possible.
    """
    x = tf.convert_to_tensor(value=x)
    with tf.compat.v1.name_scope('permute_and_reflect'):
      if self.permutable_axes.size > 0:
        x = permute_axes(x, self.full_permutation, self.permutable_axes)
      if self.reflectable_axes.size > 0:
        x = tf.reverse(x, self.reflected_axes)
      return x

max_stars_repo_name: Hanyu-Li/ffn

max_issues_repo_name: Hanyu-Li/ffn

max_forks_repo_name: Hanyu-Li/ffn


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

# 3rd party imports
import numpy as np
import xarray as xr

__author__ = "Louis Richard"
__email__ = "louisr@irfu.se"
__copyright__ = "Copyright 2020-2021"
__license__ = "MIT"
__version__ = "2.3.7"
__status__ = "Prototype"

def new_xyz(inp, trans_mat):
    r"""Transform the input field to the new frame.

    Parameters
    ----------
    inp : xarray.DataArray
        Time series of the input field in the original coordinate system.
    trans_mat : array_like
        Transformation matrix.

    Returns
    -------
    out : xarray.DataArray
        Time series of the input in the new frame.

    Examples
    --------
    >>> from pyrfu import mms, pyrf

    Time interval

    >>> tint = ["2019-09-14T07:54:00.000", "2019-09-14T08:11:00.000"]

    Spacecraft indices

    >>> mms_id = 1

    Load magnetic field and electric field

    >>> b_xyz = mms.get_data("B_gse_fgm_srvy_l2", tint, mms_id)
    >>> e_xyz = mms.get_data("E_gse_edp_fast_l2", tint, mms_id)

    Compute MVA frame

    >>> b_lmn, l, mva = pyrf.mva(b_xyz)

    Move electric field to the MVA frame

    >>> e_lmn = pyrf.new_xyz(e_xyz, mva)

    """

    if inp.data.ndim == 3:
        out_data = np.matmul(np.matmul(trans_mat.T, inp.data), trans_mat)
    else:
        out_data = (trans_mat.T @ inp.data.T).T

    out = xr.DataArray(out_data, coords=inp.coords, dims=inp.dims,
                       attrs=inp.attrs)

    return out

ext: py

max_stars_repo_name: ablotekar/irfu-python

max_issues_repo_name: ablotekar/irfu-python

max_forks_repo_name: ablotekar/irfu-python

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""
Created on Wed Jul 20 22:41:51 2016

@author: likui
"""

import csv
import numpy as np
import sys
import matplotlib.pyplot as plt #importing plotting module
from Bio.Data import IUPACData as ip
from sklearn import svm
from sklearn.metrics import roc_curve, auc

from sklearn.metrics import precision_recall_curve
from sklearn.metrics import average_precision_score
from roc import *
from sklearn import cross_validation
from sklearn.cross_validation import KFold
from sklearn.grid_search import GridSearchCV
from itertools import product
from sklearn.cross_validation import StratifiedKFold as SK
from roc import *
from numpy.random import choice
import Bio
from Bio import pairwise2
from Bio.pairwise2 import format_alignment
import numpy as np
import matplotlib.pyplot as plt #importing plotting module
from Bio import SeqIO
from temp2 import kmer, encode1,encode2
handle = open("S3.fasta", "rU")
records = list(SeqIO.parse(handle, "fasta"))
handle.close()
seq=[]
target=[]
for i in range(34):
    seq.append(str(records[i].seq))
    target.append(np.zeros((1,len(seq[i]))))

target[0][0,16:32]=1
target[0][0,87:99]=1
target[1][0,12:28]=1
target[2][0,7:20]=1
target[2][0,13:21]=1
target[2][0,19:29]=1
target[2][0,20:30]=1
target[2][0,29:37]=1
target[3][0,0:93]=1
target[4][0,59:70]=1
target[5][0,10:25]=1
target[5][0,24:35]=1
target[5][0,29:40]=1
target[5][0,36:42]=1
target[6][0,20:31]=1
target[6][0,32:41]=1
target[6][0,58:71]=1
target[6][0,82:89]=1
target[6][0,90:96]=1
target[7][0,10:20]=1
target[7][0,100:110]=1
target[7][0,115:126]=1
target[7][0,145:152]=1
target[8][0,14:20]=1
target[9][0,80:125]=1
target[10][0,0:35]=1
target[10][0,35:67]=1
target[11][0,97:103]=1
target[12][0,0:29]=1
target[12][0,100:118]=1
target[13][0,172:230]=1
target[14][0,217:289]=1
target[15][0,12:18]=1
target[16][0,10:17]=1
target[17][0,491:509]=1
target[18][0,537:545]=1
target[19][0,8:34]=1
target[20][0,4:14]=1
target[20][0,24:34]=1
target[20][0,55:61]=1
target[21][0,83:105]=1
target[21][0,104:125]=1
target[21][0,147:153]=1
target[21][0,153:163]=1
target[21][0,155:171]=1
target[21][0,179:196]=1
target[21][0,208:231]=1
target[22][0,31:41]=1
target[22][0,41:50]=1
target[23][0,65:72]=1
target[24][0,111:157]=1
target[25][0,6:21]=1
target[25][0,19:34]=1
target[25][0,42:57]=1
target[26][0,23:32]=1
target[27][0,0:142]=1
target[28][0,0:0:12]=1
target[29][0,5:12]=1
target[30][0,34:44]=1
target[30][0,48:59]=1
target[30][0,59:68]=1
target[30][0,68:82]=1
target[30][0,85:95]=1
target[31][0,588:600]=1
target[32][0,9:20]=1
target[32][0,104:115]=1
target[33][0,0:88]=1

protein=[]
labels=[]
for i in range(34):
    f=open('pro_all'+str(i)+'.txt','w')
    ww=[]#ds=np.array(3)# to store the score + labels
    w1=[]
    for yy in range(len(seq[i])-6):
        aa=seq[i][yy:yy+6]
        bb=sum(target[i][0,yy:yy+6])
       # ww=ww+[aa]
        if bb>1.0:
            w1.append(1)
            f.write(aa)
            f.write(","+str(1))
            f.write("\n")
        elif bb==0.0:
            f.write(aa)
            w1.append(0)
            f.write(","+str(0))
            f.write("\n")
        else:
            pass
    f.close()
       # f.write("\n")
      #  protein.append(ww)
       # labels.append(w1)

#test=protein[0:12]
#test_labels=labels[0:12]
#train=protein[12::]
#train_labels=labels[12::]
#score=[]
#tlb=[]
#for i in range(1):
#    test=protein[i]
#    test_labels=labels[i]
#    train=protein[0:i]+protein[i+1::]
#    train_labels=labels[0:i]+labels[i+1::]
#    Xtr=[]
#    Xla=[]
#    for ii in range(33):
#        for yy in range(len(train[ii])):
#           Xtr.append(train[ii][yy])
#           Xla.append(train_labels[ii][yy])
#
#    scr=svm.SVC(kernel='linear', class_weight='auto')
#    parameters={'C':[0.00001,0.00001,0.0001,0.001,0.01,0.1,1,10]}#,'gamma':[0.0001,0.001,0.01,0.1,1,10,100,1000]}
#    clf = GridSearchCV(scr, parameters, cv=2, scoring='roc_auc')
#    clf.fit(Xtr,Xla)
#       # ad=clf.best_params_
#    print clf.best_params_
#    print clf.best_score_
#    clff= clf.best_estimator_# initaiate the SVM instance
#    clff.fit(Xtr, Xla)
#
#
#
#    for jj in range(len(test)):
#        qq=clff.decision_function(test[jj])
#        score.append(qq)
#        tlb.append(test_labels[jj])
#
#fpr, tpr, thresholds = roc_curve(tlb,score) # plotting ROC
#a=auc(fpr,tpr)
#print a
#plt.plot(fpr,tpr,marker='.')
#plt.xlabel('fpr')
#plt.ylabel('tpr')
#plt.xlim([0,1])
#plt.ylim([0,2])
# # plt.title('ROC Linear kernel for SVM ')
#plt.grid()
#
#
#plt.legend(['auc= '+str(round(a,3))],loc=4)

reason: import numpy,from numpy


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

# This code is part of Qiskit.
#
# (C) Copyright IBM 2019.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Assertion of uniform states.
"""
from scipy.stats import chisquare
from qiskit.assertions.asserts import Asserts
from qiskit.circuit.quantumcircuit import QuantumCircuit

class AssertUniform(Asserts):
    """
        A measurement instruction that additionally performs statistical tests on the measurement
        outcomes to assert whether the state is a uniform superposition state or not.
    """
    def __init__(self, qubit, cbit, pcrit, negate):
        """
        Constructor for AssertUniform

        Args:
            qubit(QuantumRegister or list): quantum register
            cbit(ClassicalRegister or list): classical register
            pcrit(float): the critical p-value
            negate(bool): True if assertion passed is negation of statistical test passed
        """
        type_str = "Not Uniform" if negate else "Uniform"
        super().__init__(self._syntax_for_measure(qubit), self._syntax_for_measure(cbit),
                         pcrit, negate, type_str)

    def stat_test(self, counts):
        """
        Performs a chi-squared statistical test on the experimental outcomes.  Internally, compares
        a normalized table of experimental counts to the scipy.stats.chisquare default, for which
        all outcomes are equally likely.

        Args:
            counts(dictionary): result.get_counts(experiment)

        Returns:
            tuple: tuple containing:

                chisq(float): the chi-square value

                pval(float): the p-value

                passed(Boolean): if the test passed
        """

        vals_list = list(counts.values())

        # account for bitstring with zero counts
        numqubits = len(list(counts)[0])
        numzeros = 2 ** numqubits - len(counts)
        vals_list.extend([0] * numzeros)

        chisq, pval = chisquare(vals_list)
        passed = bool(pval >= self._pcrit)
        return (chisq, pval, passed)

def get_breakpoint_uniform(self, qubit, cbit, pcrit=0.05):
    """
    Creates a breakpoint, which is a renamed deep copy of the QuantumCircuit, and creates and
    appends an AssertUniform instruction to its end.  If the statistical test passes, the
    assertion passes; if the test fails, the assertion fails.

    Args:
        qubit (QuantumRegister or list): quantum register
        cbit (ClassicalRegister or list): classical register
        pcrit (float): critical p-value for the hypothesis test

    Returns:
        QuantumCircuit: copy of quantum circuit at the assert point
    """
    clone = self.copy(Asserts._new_breakpoint_name())
    assertion = AssertUniform(qubit, cbit, pcrit, False)
    clone.append(assertion, [assertion._qubit], [assertion._cbit])
    return clone

QuantumCircuit.get_breakpoint_uniform = get_breakpoint_uniform

def get_breakpoint_not_uniform(self, qubit, cbit, pcrit=0.05):
    """
    Creates a breakpoint, which is a renamed deep copy of the QuantumCircuit, and creates and
    appends an AssertUniform instruction to its end.  If the statistical test passes, the
    assertion fails; if the test fails, the assertion passes.

    Args:
        qubit (QuantumRegister or list): quantum register
        cbit (ClassicalRegister or list): classical register
        pcrit (float): critical p-value for the hypothesis test

    Returns:
        QuantumCircuit: copy of quantum circuit at the assert point
    """
    clone = self.copy(Asserts._new_breakpoint_name())
    assertion = AssertUniform(qubit, cbit, pcrit, True)
    clone.append(assertion, [assertion._qubit], [assertion._cbit])
    return clone

QuantumCircuit.get_breakpoint_not_uniform = get_breakpoint_not_uniform

max_stars_repo_name: edasgupta/qiskit-terra

max_issues_repo_name: edasgupta/qiskit-terra

max_forks_repo_name: edasgupta/qiskit-terra

reason: from scipy


text:
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import norm
from lightgbm import LGBMClassifier
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from  sklearn.preprocessing import RobustScaler, StandardScaler, Normalizer, QuantileTransformer, MinMaxScaler
from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier
import optuna
from sklearn.metrics import accuracy_score

import os
for dirname, _, filenames in os.walk('/kaggle/input'):
    for filename in filenames:
        print(os.path.join(dirname, filename))
train=pd.read_csv('/kaggle/input/forest-cover-type-prediction/train.csv', index_col='Id')
test=pd.read_csv('/kaggle/input/forest-cover-type-prediction/test.csv', index_col='Id')
submission=pd.read_csv('/kaggle/input/forest-cover-type-prediction/sampleSubmission.csv', index_col='Id')
Ytrain=train['Cover_Type']
train=train[list(test)]
train.isnull().sum()
test.isnull().sum()
for col in ['Wilderness_Area1', 'Wilderness_Area2', 'Wilderness_Area3', 'Wilderness_Area4', 'Soil_Type1', 'Soil_Type2', 'Soil_Type3', 'Soil_Type4', 'Soil_Type5', 'Soil_Type6',
            'Soil_Type7', 'Soil_Type8', 'Soil_Type9', 'Soil_Type10', 'Soil_Type11', 'Soil_Type12', 'Soil_Type13', 'Soil_Type14', 'Soil_Type15', 'Soil_Type16', 'Soil_Type17',
            'Soil_Type18', 'Soil_Type19', 'Soil_Type20', 'Soil_Type21', 'Soil_Type22', 'Soil_Type23', 'Soil_Type24', 'Soil_Type25', 'Soil_Type26', 'Soil_Type27', 'Soil_Type28',
            'Soil_Type29', 'Soil_Type30', 'Soil_Type31', 'Soil_Type32', 'Soil_Type33', 'Soil_Type34', 'Soil_Type35', 'Soil_Type36', 'Soil_Type37', 'Soil_Type38', 'Soil_Type39',
            'Soil_Type40']:
    if train[col].nunique() != test[col].nunique():
        print(f'Different Values train-test : {col}')
train=train.drop(columns=['Soil_Type7','Soil_Type15'])
test=test.drop(columns=['Soil_Type7','Soil_Type15'])
sns.countplot(Ytrain)
Ytrain.value_counts()
train.head()
train['Distance_To_Hydrology']=np.sqrt((train['Horizontal_Distance_To_Hydrology'] **2)  + (train['Vertical_Distance_To_Hydrology'] **2))
test['Distance_To_Hydrology']=np.sqrt((test['Horizontal_Distance_To_Hydrology'] **2)  + (test['Vertical_Distance_To_Hydrology'] **2))

figure, (ax1,  ax2) = plt.subplots(nrows=1, ncols=2)
figure.set_size_inches(14,6)

sns.pointplot(data=train, x=Ytrain, y='Distance_To_Hydrology', ax=ax1)
sns.distplot(train['Distance_To_Hydrology'], ax=ax2)
print(f"Num_Kind of plants living in water? : {Ytrain.loc[train.loc[train['Distance_To_Hydrology']==0].index].nunique()}")
sns.distplot(np.log(train['Distance_To_Hydrology']+1), fit=norm)
train['Distance_To_Hydrology']=np.log(train['Distance_To_Hydrology']+1)
test['Distance_To_Hydrology']=np.log(test['Distance_To_Hydrology']+1)

train.loc[train['Distance_To_Hydrology']==0,'Distance_To_Hydrology']=np.nan
test.loc[test['Distance_To_Hydrology']==0,'Distance_To_Hydrology']=np.nan
qut=QuantileTransformer(output_distribution='normal')
pipeline=make_pipeline(qut, PCA(n_components=1))
train.corr()
sns.scatterplot(data=train, x='Elevation', y='Horizontal_Distance_To_Roadways', hue=Ytrain)
train['Elevation_Roadways']=pipeline.fit_transform(train[['Elevation','Horizontal_Distance_To_Roadways']])
test['Elevation_Roadways']=pipeline.transform(test[['Elevation','Horizontal_Distance_To_Roadways']])
alpha=train[['Hillshade_9am','Hillshade_3pm']]
alpha2=test[['Hillshade_9am','Hillshade_3pm']]
train['Hillshade9am_Hillshade_3pm']=pipeline.fit_transform(alpha)
test['Hillshade9am_Hillshade_3pm']=pipeline.transform(alpha2)
sns.lineplot(data=train, x=Ytrain, y='Hillshade9am_Hillshade_3pm')
alpha=train[['Aspect','Hillshade9am_Hillshade_3pm']]
alpha2=test[['Aspect','Hillshade9am_Hillshade_3pm']]
train['Aspect_Hillshade']=pipeline.fit_transform(alpha)
test['Aspect_Hillshade']=pipeline.transform(alpha2)
alpha=train[['Slope','Hillshade_Noon']]
alpha2=test[['Slope','Hillshade_Noon']]
train['Slope_Hillshade']=pipeline.fit_transform(alpha)
test['Slope_Hillshade']=pipeline.transform(alpha2)
Xtrain=train
Xtest=test
rf=RandomForestClassifier(n_estimators=5000, random_state=18, n_jobs=-1)
rf.fit(Xtrain.fillna(-999), Ytrain)
feature_importance_df=pd.DataFrame(data=None,  columns=['feature','importances'])
feature_importance_df['importances']=rf.feature_importances_
feature_importance_df['feature']=Xtrain.columns

feature_importance_df=feature_importance_df.sort_values(by='importances', ascending=False)
feature_importance_df.tail()
feature_names=feature_importance_df.head(56)['feature']
Xtrain=train[feature_names]
Xtest=test[feature_names]
print(Xtrain.shape, Ytrain.shape, Xtest.shape)
params={'n_estimators': 1222, 'learning_rate': 0.07307234151834806, 'num_leaves': 96, 'colsample_bytree': 0.8972376156262298,
        'subsample': 0.9312856106293543, 'min_child_samples': 1}

lightgbm=LGBMClassifier(random_state=18, subsample_freq=1, silent=False, **params)
lightgbm.fit(Xtrain, Ytrain)
predictions=lightgbm.predict(Xtest)
submission['Cover_Type']=predictions
submission.to_csv('LGBSingleModel.csv')
submission.head()

max_stars_repo_name: josepablocam/janus-public

max_issues_repo_name: josepablocam/janus-public

max_forks_repo_name: josepablocam/janus-public

reason: import numpy,from scipy


text:
import pandas as pd
import numpy as np
import requests

# class name，必須跟檔案名一致，例如 class INSTITUTIONAL_INVESTORS_FOLLOWER，檔名也是 INSTITUTIONAL_INVESTORS_FOLLOWER.py
class INSTITUTIONAL_INVESTORS_FOLLOWER:
    def __init__(self,
                 stock_price,
                 **kwargs, ):
        # -------------------------------------------------------------------
        # 此區塊請勿更動
        stock_price = stock_price.sort_values('date')
        # 股價
        self.stock_price = stock_price
        # 融資融券
        self.MarginPurchaseShortSale = kwargs.get("MarginPurchaseShortSale", pd.DataFrame())
        # 三大法人買賣
        self.InstitutionalInvestorsBuySell = kwargs.get("InstitutionalInvestorsBuySell", pd.DataFrame())
        # 外資持股
        self.Shareholding = kwargs.get("Shareholding", pd.DataFrame())
        # 此區塊請勿更動
        # -------------------------------------------------------------------
        self.url = 'https://www.finlab.tw/%E8%85%A6%E5%8A%9B%E6%BF%80%E7%9B%AA%E7%9A%84%E5%A4%96%E8%B3%87%E7%AD%96%E7%95%A5%EF%BC%81/'
        self.summary = '''策略概念: 法人大量買超會導致股價上漲, 賣超反之
                          策略規則: 三大法人大量買超隔天就賣，大量賣超就買'''
        self.InstitutionalInvestorsBuySell[['sell', 'buy']] = self.InstitutionalInvestorsBuySell[
            ['sell', 'buy']].astype(int)
        self.InstitutionalInvestorsBuySell = self.InstitutionalInvestorsBuySell.groupby(['date', 'stock_id'],
                                                                                as_index=False).agg(
            {'buy': np.sum, 'sell': np.sum})
        self.InstitutionalInvestorsBuySell['diff'] = self.InstitutionalInvestorsBuySell['buy'] - \
                                                     self.InstitutionalInvestorsBuySell['sell']

        self.stock_price = pd.merge(self.stock_price[['stock_id', 'date', 'open', 'close']],
                                    self.InstitutionalInvestorsBuySell[['stock_id', 'date', 'diff']],
                                    on=['stock_id', 'date'], how='left').fillna(0)
        self.stock_price['signal'] = self.detect_Abnormal_Peak(
            y=self.stock_price['diff'].values, lag=10, threshold=3, influence=0.35)['signals']

    def detect_Abnormal_Peak(self, y, lag, threshold, influence):
        signals = np.zeros(len(y))
        filteredY = np.array(y)
        avgFilter = [0] * len(y)
        stdFilter = [0] * len(y)
        avgFilter[lag - 1] = np.mean(y[0:lag])
        stdFilter[lag - 1] = np.std(y[0:lag])
        for i in range(lag, len(y)):
            if abs(y[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1]:
                if y[i] > avgFilter[i - 1]:
                    signals[i] = 1
                else:
                    signals[i] = -1

                filteredY[i] = influence * y[i] + (1 - influence) * filteredY[i - 1]
                avgFilter[i] = np.mean(filteredY[(i - lag + 1):i + 1])
                stdFilter[i] = np.std(filteredY[(i - lag + 1):i + 1])
            else:
                signals[i] = 0
                filteredY[i] = y[i]
                avgFilter[i] = np.mean(filteredY[(i - lag + 1):i + 1])
                stdFilter[i] = np.std(filteredY[(i - lag + 1):i + 1])

        return dict(signals=np.asarray(signals),
                    avgFilter=np.asarray(avgFilter),
                    stdFilter=np.asarray(stdFilter))

    def trade(self, date):
        '''
        此區塊，可進行資料處理、做技術指標，寫自己的策略，
        寫你自己的策略, 必須 return : 1 (買) or -1 (賣) or 0 (不操作)
        根據時間date，回傳當下要進行什麼操作 ( 買/賣/不操作 )

        date : 昨天時間
        用昨天的資料，計算技術指標，判斷今天買/賣
        '''
        # example
        value = self.stock_price[self.stock_price['date'] == date]
        signal = value['signal'].values
        if len(value) == 0:
            return 0

        if signal == -1:
            return 1  # buy
        elif signal == 1:
            return -1  # sell
        elif signal == 0:
            return 0  # no buy and no sell

def test():
    stock_id = '0056'
    date = '2015-01-01'

    url = 'http://finmindapi.servebeer.com/api/data'
    form_data = {'dataset': 'TaiwanStockPrice',
                 'stock_id': stock_id,
                 'date': date}

    res = requests.post(url, verify=True, data=form_data)

    temp = res.json()
    stock_price = pd.DataFrame(temp['data'])

    form_data = {'dataset': 'InstitutionalInvestorsBuySell',
                 'stock_id': stock_id,
                 'date': date}
    res = requests.post(
        url, verify=True,
        data=form_data)

    temp = res.json()
    InstitutionalInvestorsBuySell = pd.DataFrame(temp['data'])

    obj = INSTITUTIONAL_INVESTORS_FOLLOWER(
        stock_price=stock_price,
        InstitutionalInvestorsBuySell=InstitutionalInvestorsBuySell
    )
    obj.trade('2019-05-07')

max_stars_repo_name: Jacarlianda/FinMind

max_issues_repo_name: Jacarlianda/FinMind

max_forks_repo_name: Jacarlianda/FinMind


text:
#%%
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import phd.viz
colors, palette = phd.viz.phd_style()

# Load in the datasets
data_a = pd.read_csv('../../data/ch6_induction_si/figS1_partA.csv')
data_b = pd.read_csv('../../data/ch6_induction_si/figs1_partB.csv')
data_O2 = pd.read_csv('../../data/ch2_induction/RazoMejia_2018.csv', comment='#')
data_O2['repressors'] *= 2
data_O2 = data_O2[(data_O2['operator'] == 'O2') & (data_O2['repressors'] == 260)]

# Define necessary parameters.
ka = data_a[data_a.parameter == 'logKA']
ki = data_a[data_a.parameter == 'logKI']
c_range = np.logspace(-8, -2, 500)
R = 260
ep_r = -13.9
_colors = sns.color_palette('magma', n_colors=len(data_b)+1)

def foldchange(c, R, ep_ai, ep_r, ep_a, ep_i):
    mwc_term = (1 + c / ep_a)**2 * (1 + np.exp(-ep_ai)) / ((1 + c / ep_a)**2 +
                                                           np.exp(-ep_ai) * (1 + c / ep_i)**2)
    fc = (1 + mwc_term * (R / 4.6E6) * np.exp(-ep_r))**-1
    return fc
# Define the figure axis and labels.
fig, ax = plt.subplots(1, 2, figsize=(6, 2.25))
phd.viz.despine(ax)
_ = ax[0].set_xlabel(r'allosteric parameter $\Delta\varepsilon_{AI}\,(k_BT)$')
_ = ax[0].set_ylabel(r'best-fit parameter value')
_ = ax[1].set_xscale('log')
_ = ax[1].set_xlabel('IPTG [µM]')
_ = ax[1].set_ylabel('fold-change')

# Plots for panel (A)
_ = ax[0].plot(ka.ep, ka.bestfit, '-', color=colors['orange'],
               label=r'$\mathrm{log}\, \frac{K_A}{1\mathrm{M}}$')
_ = ax[0].plot(ki.ep, ki.bestfit, '-', color=colors['purple'],
               label=r'$\mathrm{log}\, \frac{K_A}{1\mathrm{M}}$')

# Plot the curves
for i in range(len(data_b)):
    ep_ai = data_b.iloc[i]['ep_ai']
    ka = np.exp(data_b.iloc[i]['log_ka'])
    ki = np.exp(data_b.iloc[i]['log_ki'])
    fc = foldchange(c_range, R, ep_ai, ep_r, ka, ki)
    _ = ax[1].plot(c_range * 1E6, fc, label=ep_ai, color=_colors[i])

# plot the data.
grouped = data_O2.groupby(['IPTG_uM']).fold_change_A
for group, data in grouped:
    mean_fc = np.mean(data)
    mean_sem = np.std(data) / np.sqrt(len(data))
    _ = ax[1].errorbar(group, mean_fc, mean_sem,
                       linestyle='none', color=colors['orange'],
                       fmt='o', markeredgecolor='white', markeredgewidth=0.5,
                       label='__nolegend__', ms=4.5)

# Add the legends and labels.
_ = ax[0].legend(loc='lower left')
leg = ax[1].legend(loc='upper left', title=r"""allosteric parameter
       $\Delta\varepsilon_{AI}$ $(k_BT)$""", bbox_to_anchor=(1, 1),
       fontsize=6)
leg.get_title().set_fontsize(6)
fig.text(0, 0.95, '(A)', fontsize=8)
fig.text(0.45, 0.95, '(B)', fontsize=8)

# Format and save the figure.
plt.tight_layout()
plt.savefig('../figs/ch6_figS1.pdf', bbox_inches='tight')
plt.savefig('../figs/ch6_figS1.png', bbox_inches='tight')

# %%

max_stars_repo_name: gchure/phd

max_issues_repo_name: gchure/phd

max_forks_repo_name: gchure/phd


text:
import numpy as np
from scipy import sparse, stats
import sklearn.utils.sparsefuncs as sf

def diffexp_ttest(meanA,vA,nA,meanB,vB,nB, top_n=8, diffexp_lfc_cutoff=0.01):
    return diffexp_ttest_from_mean_var(meanA, vA, nA, meanB, vB, nB, 1000, diffexp_lfc_cutoff)

def diffexp_ttest_from_mean_var(meanA, varA, nA, meanB, varB, nB, top_n, diffexp_lfc_cutoff):
    n_var = meanA.shape[0]
    top_n = min(n_var,top_n)

    # variance / N
    vnA = varA / min(nA, nB)  # overestimate variance, would normally be nA
    vnB = varB / min(nA, nB)  # overestimate variance, would normally be nB
    sum_vn = vnA + vnB

    # degrees of freedom for Welch's t-test
    with np.errstate(divide="ignore", invalid="ignore"):
        dof = sum_vn ** 2 / (vnA ** 2 / (nA - 1) + vnB ** 2 / (nB - 1))
    dof[np.isnan(dof)] = 1

    # Welch's t-test score calculation
    with np.errstate(divide="ignore", invalid="ignore"):
        tscores = (meanA - meanB) / np.sqrt(sum_vn)
    tscores[np.isnan(tscores)] = 0

    # p-value
    pvals = stats.t.sf(np.abs(tscores), dof) * 2
    pvals_adj = pvals * n_var
    pvals_adj[pvals_adj > 1] = 1  # cap adjusted p-value at 1

    # logfoldchanges: log2(meanA / meanB)
    logfoldchanges = np.log2(np.abs((meanA + 1e-9) / (meanB + 1e-9)))

    stats_to_sort = tscores
    # find all with lfc > cutoff
    lfc_above_cutoff_idx = np.nonzero(np.abs(logfoldchanges) > diffexp_lfc_cutoff)[0]

    # derive sort order
    sort_order = np.argsort(stats_to_sort)
    sort_order = np.concatenate((sort_order[-top_n:][::-1],sort_order[:top_n][::-1]))
    if lfc_above_cutoff_idx.shape[0] > 0:
        sort_order = sort_order[np.in1d(sort_order,lfc_above_cutoff_idx)]
    # top n slice based upon sort order
    logfoldchanges_top_n = logfoldchanges[sort_order]
    pvals_top_n = pvals[sort_order]
    pvals_adj_top_n = pvals_adj[sort_order]

    # varIndex, logfoldchange, pval, pval_adj
    result = {"positive": [[sort_order[i], logfoldchanges_top_n[i], pvals_top_n[i], pvals_adj_top_n[i]] for i in
                           range(top_n)],
              "negative": [[sort_order[i], logfoldchanges_top_n[i], pvals_top_n[i], pvals_adj_top_n[i]] for i in
                           range(-1, -1 - top_n, -1)], }

    return result

def mean_var_n(X):
    """
    Two-pass variance calculation.  Numerically (more) stable
    than naive methods (and same method used by numpy.var())
    https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Two-pass
    """
    n = X.shape[0]
    if sparse.issparse(X):
        mean,v = sf.mean_variance_axis(X,axis=0)
    else:
        mean,v = X.mean(0),X.var(0)
    meansq = v - mean**2
    return mean, meansq, v, n

max_stars_repo_name: czbiohub/cellxgene

max_issues_repo_name: czbiohub/cellxgene

max_forks_repo_name: czbiohub/cellxgene


text:
import numpy as np
import hashlib
import json
from datetime import datetime
from shapely.geometry import Point

def translate_coordinate(lat, lng, l, h):
    """
    This function takes in a latitude,longitude pair and translates it to produce a l x h area (m^2) with the
    original latitude, longitude pair located at the centroid.
    """
    # Define the Earth’s radius, assuming a spherical surface
    earth_radius = 6378137.0

    # Translate location in meters
    dn = h / 2.0
    de = l / 2.0

    # Coordinate translation in radians
    dLat = dn / earth_radius
    dLng = de / (earth_radius * np.cos(np.pi * lat / 180.0))

    # Translate position in decimal degrees
    max_lat = lat + dLat * 180.0 / np.pi
    max_lng = lng + dLng * 180.0 / np.pi
    min_lat = lat - dLat * 180.0 / np.pi
    min_lng = lng - dLng * 180.0 / np.pi

    return max_lat, max_lng, min_lat, min_lng

def identify_centroid(max_lat, max_lng, min_lat, min_lng):
    """
    This function takes in the coordinates of the bounding box and outputs its centroids.
    """
    return (max_lat + min_lat) / 2.0, (max_lng + min_lng) / 2.0

def calculate_circle_radius(max_lat, max_lng, centre_lat, centre_lng):
    """
    This function takes in the edges of a bounding box and its centroid to calculate the radius of a bounding
    circle that can be fitted within the bounding box with the centroid at its centre.
    """
    max_coordinate_dist_lat = abs(max_lat - centre_lat)
    max_coordinate_dist_lng = abs(max_lng - centre_lng)

    # Define the Earth’s radius, assuming a spherical surface
    earth_radius = 6378137.0

    lng_radius = max_coordinate_dist_lng * (earth_radius * np.cos(np.pi * centre_lat / 180.0)) / (180.0 / np.pi)
    lat_radius = max_coordinate_dist_lat * earth_radius / (180.0 / np.pi)

    if lng_radius >= lat_radius:
        print('Radius of bounding circle (Latitude): {}'.format(lat_radius))
        return lat_radius
    else:
        print('Radius of bounding circle (Longitude): {}'.format(lng_radius))
        return lng_radius

def calculate_lat_lng_distance(min_coordinate, max_coordinate, is_lat=True, latitude=None):
    """
    This function calculates the distance (in metres) between two coordinates. Calculations differ depending if the
    coordinate pair comes from the latitude or longitude.
    """
    # Define the Earth’s radius, assuming a spherical surface
    radius = 6378137.0

    if is_lat:
        return (max_coordinate - min_coordinate) * radius / (180.0 / np.pi)
    else:
        return (max_coordinate - min_coordinate) * (radius * np.cos(np.pi * latitude / 180.0)) / (180.0 / np.pi)

def generate_coordinate_list(min_coordinate, max_coordinate, num_box):
    """
    This function takes in a pair of coordinates (either from the latitude or longitude) and breaks the entire range
    into multiple intervals, whereby each interval is the size of one side of the query box. The coordinates marking
    each interval is then stored in a list and returned as an output.
    """
    query_range = max_coordinate - min_coordinate
    query_interval_length = query_range / num_box

    coordinate_list = []
    coordinate = min_coordinate
    for i in range(int(num_box) + 1):
        coordinate_list.append(coordinate)
        coordinate += query_interval_length

    return coordinate_list

def generate_bounding_coordinates(lat_list, lng_list):
    """
    This function takes in two lists containing coordinates in the latitude and longitude direction and generates a
    list containing coordinates for each bounding box.
    """
    lat_list = list(reversed(lat_list))
    coordinate_list = []
    for i in range(len(lng_list) - 1):
        for j in range(len(lat_list) - 1):
            coordinate_list.append([lat_list[j + 1], lng_list[i], lat_list[j], lng_list[i + 1]])

    return coordinate_list

def divide_bounding_box(max_lat, min_lat, max_lng, min_lng, querybox_dim):
    """
    This function divides the bounding box (defined using its coordinates) up into smaller sub-bounding boxes and
    returns a list of query box coordinates.
    """
    lat_dist = calculate_lat_lng_distance(min_coordinate=min_lat, max_coordinate=max_lat)
    lng_dist = calculate_lat_lng_distance(min_coordinate=min_lng, max_coordinate=max_lng, is_lat=False,
                                          latitude=(max_lat + min_lat) / 2.0)
    num_box_lat = round(lat_dist / querybox_dim)
    num_box_lng = round(lng_dist / querybox_dim)
    lat_list = generate_coordinate_list(min_coordinate=min_lat, max_coordinate=max_lat, num_box=num_box_lat)
    lng_list = generate_coordinate_list(min_coordinate=min_lng, max_coordinate=max_lng, num_box=num_box_lng)
    lat_lng_list = generate_bounding_coordinates(lat_list=lat_list, lng_list=lng_list)

    return lat_lng_list

def pixelise_region(coordinates, shapefile):
    """
    This function filters out a list of coordinates based on whether it intersects with the regions stored within
    the shapefile.
    """
    return [coordinate for coordinate in coordinates if
            (np.sum(shapefile['geometry'].apply(lambda x: Point(coordinate[1], coordinate[0]).within(x))) != 0) |
            (np.sum(shapefile['geometry'].apply(lambda x: Point(coordinate[3], coordinate[0]).within(x))) != 0) |
            (np.sum(shapefile['geometry'].apply(lambda x: Point(coordinate[1], coordinate[2]).within(x))) != 0) |
            (np.sum(shapefile['geometry'].apply(lambda x: Point(coordinate[3], coordinate[2]).within(x))) != 0)]

def generate_id(poi_dict):
    """
    Generates an unique ID based on the encoded string of the POI's name, address and location information.
    """
    poi_string = poi_dict['properties']['name'] + \
                 poi_dict['properties']['address'] + \
                 str(poi_dict['geometry']['lat']) + \
                 str(poi_dict['geometry']['lng'])
    hasher = hashlib.md5()
    hasher.update(poi_string.encode('utf-8'))
    generated_id = hasher.hexdigest()
    return generated_id

def within_boundary_area(lat, lng, min_lat, max_lat, min_lng, max_lng):
    """
    Checks if the latitude and longitude pair falls within a region of interest.
    """
    if lat > max_lat or lat < min_lat:
        return False
    if lng > max_lng or lng < min_lng:
        return False

    return True

def extract_date():
    """
    Extract the date when the query is made in the format YYYYMMDD.
    """
    return datetime.today().strftime('%Y%m%d')

def remove_duplicate(filename):
    """
    Remove duplicated POIs in the dataset using the POI's unique ID.
    """
    # load json file from local directory
    with open(filename) as json_file:
        feature_collection = json.load(json_file)

    # initialise storage for duplicated IDs
    dropped_index = []
    id_list = []

    for i in range(len(feature_collection['features'])):
        if feature_collection['features'][i]['id'] in id_list:
            dropped_index.append(i)
        else:
            id_list.append(feature_collection['features'][i]['id'])

    for index in sorted(dropped_index, reverse=True):
        del feature_collection['features'][index]

    # save json file on local directory
    with open(filename, 'w') as json_file:
        json.dump(feature_collection, json_file)

def capitalise_string(string):
    """
    Capitalise the first letter of each word in a string. The original string may contain ().
    """
    capitalised_string = ''
    string_list = string.lower().split(' ')
    for i in range(len(string_list)):
        if not string_list[i]:
            continue
        elif string_list[i][0] != '(':
            capitalised_string += string_list[i].capitalize() + ' '
        else:
            capitalised_string += '(' + string_list[i][1:].capitalize() + ' '

    return capitalised_string[:-1]


text:
# SPDX-License-Identifier: Apache-2.0

import logging
from fed_learn.client.fed_privacy import PrivacyProtocol
import numpy as np

class MyPrivacyProtocol(PrivacyProtocol):

    def __init__(self, percentile=10, gamma=0.01):
        self.logger = logging.getLogger('---- my Pricavy protocol')

        PrivacyProtocol.__init__(self)

    def apply(self, model_diff, train_ctx=None):
        self.logger.info("----------------")
        self.logger.info("-------client applying privacy---------")
        self.logger.info("----------------")
        """
        :param model_diff: model after a round of local training
        :param train_ctx: unused param.
        :return: model data to be uploaded to the server
        """

        # invariant to local steps
        scale = np.float(train_ctx.total_steps)
        delta_w = {name: model_diff[name] / scale for name in model_diff}

        return delta_w

class MyRandomProtocol(PrivacyProtocol):
    """
    Randomly drop some gradients
    """

    def __init__(self, percentage=10, gamma=0.01, non_zero_only=True):
        self.logger = logging.getLogger('RandomProtocol')

        PrivacyProtocol.__init__(self)
        # must be in 0..100, only update abs diff greater than percentile
        self.percentage = percentage
        assert percentage >= 0 & percentage <= 100, "select a percentage between 0 and 100!"

    def apply(self, model_diff, train_ctx=None):
        delta_w = {}
        n_removed = 0
        n_total = 0
        for name in model_diff:
            diff = model_diff[name]
            n = np.size(diff)
            n_total += n
            n_select = int(np.floor(n*self.percentage/100))
            if n_select > 0:
                idx = np.arange(n)
                select = np.random.choice(idx, size=n_select, replace=False)
                np.put(diff, select, v=0.0)  # set the randomly selected gradients to zero
                n_removed += n_select
            delta_w[name] = diff

        self.logger.info(f'Removed {n_removed} of {n_total} ({100*n_removed/n_total:.2f}%) parameters.')
        return delta_w

max_stars_repo_name: ZQyou/clara-train-examples

max_issues_repo_name: ZQyou/clara-train-examples

max_forks_repo_name: ZQyou/clara-train-examples


text:
import unittest

import numpy as np
from numpy.testing import assert_almost_equal, assert_array_almost_equal

from structured import common, gibbs, objects, varapprox
from structured.common import mult_diag

def _random_model(T, dy, dc, dsc, params='noninformative', w_real=None):
    # params: 'noninformative', 'exact',
    #          dictionary with tT, p00, p11, tlbd, lbds

    # w_mean set to the correct weights
    if w_real is None:
        w_real = common.random_weights_orth(dy, dc, dsc)
    T_real = np.array([[0.8,0.2],[0.1,0.9]], dtype='d')
    lbd_real = np.linspace(0.8, 0.2, dsc)

    if params=='noninformative':
        tT, p00, p11 = 1., 0.5, 0.5
        tlbd, lbds = 1., np.linspace(0.0,0.0,dsc)
        tnoise, noise_prec = 1., 1./0.01
        gamma = 1.
    elif params=='exact':
        tT, p00, p11 = 1e3, T_real[0,0], T_real[1,1]
        tlbd, lbds = 1e3, lbd_real.copy()
        tnoise, noise_prec = 1e5, 1e10
        gamma = 1.e-5

    w_mean = w_real.copy()
    w_var = np.zeros((dy,dc,dsc), dtype='d') + 1e-10
    params = objects.Parameters(dy, dc, dsc, 0.5, tT, p00, p11,
                                tlbd, lbds, tnoise, noise_prec,
                                gamma=gamma,
                                w_mean=w_mean, w_var=w_var,
                                w_real=w_real,
                                T_real = T_real, lbd_real=lbd_real)
    y, c, s = params.sample_all(T, noise=False, s_off=True)
    state = objects.State(y, params)
    state.set_to_prior()

    return params, state, c, s

class TestVarapprox(unittest.TestCase):
    def test_e_step(self):
        params, state, c, s = _random_model(100, 18, 5, 3, params='exact')

        varapprox.e_step(state, params)

        assert_array_almost_equal(c, state.pc, 7)
        # s_mn is zero if content off
        for i in range(params.dc):
            state.s_mn[:,i,:] *= np.repeat(c[:,i,np.newaxis], (params.dsc,), 1)
        assert_array_almost_equal(s, state.s_mn, 7)
        return params, state, c, s

    def test_m_step(self):
        T, dy, dc, dsc = 1000, 14, 4, 3
        params, state, c, s = _random_model(T, dy, dc, dsc, params='exact')

        params.w_mean = common.random_weights_orth(dy, dc, dsc)
        params.w_var = np.zeros((dy,dc,dsc), dtype='d') + 1e-1
        params.init_wDw()

        # vars from the data
        state.set_to(c,s)
        varapprox.close_e_step(params, state)

        # learn w
        varapprox.m_step(state, params, max_iter=1000)
        # evaluate
        assert_array_almost_equal(params.w_real, params.w_mean, 8)

    def test_m_step_rho(self):
        T, dy, dc, dsc = 5000, 14, 3, 2
        params, state, c, s = _random_model(T, dy, dc, dsc, params='exact')

        tnoise, noise_prec = 1., 25.
        params.init_noise_prior(tnoise, noise_prec)

        # add structured noise to the observations
        noise_std = np.linspace(1., 0.01, dy)*0.2*state.y.std()
        state.y += mult_diag(noise_std, np.random.normal(size=(T, dy)), left=False)
        # vars from the data
        state.set_to(c,s)
        varapprox.close_e_step(params, state)

        # learn w
        varapprox.m_step(state, params, max_iter=1000)
        # evaluate
        assert_array_almost_equal(np.sqrt(1./params.rho_mn), noise_std, 2)
        return params, state, c, s

    def test_perturbation(self):
        T, dy, dc, dsc = 2000, 14, 3, 2
        params, state, c, s = _random_model(T, dy, dc, dsc)

        params.w_mean += np.random.normal(scale=1e-1, size=(dy,dc,dsc))
        params.w_var = np.zeros((dy,dc,dsc), dtype='d') + 1e-3
        params.init_wDw()

        old_free = -np.inf
        for i in range(10):
            print 'iteration %d/10' % i
            for j in range(3):
                varapprox.e_step(state, params)
                varapprox.close_e_step(params, state)
                new_free = varapprox.free_energy(params, state)
                self.assertGreaterEqual(new_free-old_free, 0.0)
                old_free = new_free

            varapprox.close_e_step(params, state)
            varapprox.m_step(state, params, max_iter=2000)
            new_free = varapprox.free_energy(params, state)
            self.assertGreaterEqual(new_free-old_free, 0.0)

            vols, dimrel, bestc = common.evaluate_style_subspace_perms(params.w_real, params.w_mean)
            print vols

        # evaluate
        vols, dimrel, bestc = common.evaluate_style_subspace_perms(params.w_real, params.w_mean)
        self.assertLess(max(vols), 1e-2)
        assert_array_almost_equal(params.lbd_mn.mean(axis=0), params.lbd_real[0,:], 1)
        assert_array_almost_equal(params.T_mn, params.T_real, 1)

    def test_learn_gamma(self):
        dy, dc, dsc = 7, 2, 3
        w_real = common.random_weights_orth(dy, dc, dsc) * np.sqrt(dy)
        w_real[:,0,0] = 0.
        w_real[:,1,0] *= np.sqrt(2)
        params, state, c, s = _random_model(1000, dy, dc, dsc, params='exact',
                                            w_real=w_real)

        # vars from the data
        state.set_to(c,s)

        varapprox.learn_gamma(params)
        self.assertGreater(params.gamma[0,0], 1e4)
        assert_almost_equal(params.gamma[1,0], 0.5, 8)
        for i in range(dc):
            for j in range(1,dsc):
                assert_almost_equal(params.gamma[i,j], 1., 8)

    def test_sample_c(self):
        params, state, c, s = _random_model(50, 10, 2, 3, params='exact')
        y, c, s = params.sample_all(50, noise=False, s_off=False)
        state = objects.State(y, params)
        state.set_to_prior()

        nsamples = 300
        samples = np.zeros((nsamples+1, 50, 2))
        samples[0,:] = state.pc.copy()
        for i in range(nsamples):
            samples[i+1,:] = gibbs.sample_c(state.y, samples[i,:], s, params)

        assert_almost_equal(c, samples.mean(axis=0), 2)

    def test_sample_s(self):
        T, dy, dc, dsc = 50, 10, 3 ,2
        params, state, c, s = _random_model(T, dy, dc, dsc, params='noninformative')

        nsamples = 300
        s_samples = np.zeros((nsamples+1, T, dc, dsc))
        s_samples[0,:] = state.s_mn.copy()
        for i in range(nsamples):
            s_samples[i+1,:] = gibbs.sample_s(state.y, c, s_samples[i,:], params)

        assert_almost_equal(s, s_samples.mean(axis=0), 1)

max_stars_repo_name: pberkes/identity-attributes-model

max_issues_repo_name: pberkes/identity-attributes-model

max_forks_repo_name: pberkes/identity-attributes-model


text:
from hoomd.snapshot import Snapshot
import hoomd
import numpy
import pytest

@pytest.fixture(scope='function')
def snap(device):
    s = Snapshot(device.communicator)
    N = 1000

    if s.communicator.rank == 0:
        s.configuration.box = [20, 20, 20, 0, 0, 0]

        s.particles.N = N
        s.particles.position[:] = numpy.random.uniform(-10, 10, size=(N, 3))
        s.particles.velocity[:] = numpy.random.uniform(-1, 1, size=(N, 3))
        s.particles.typeid[:] = numpy.random.randint(0, 3, size=N)
        s.particles.mass[:] = numpy.random.uniform(1, 2, size=N)
        s.particles.charge[:] = numpy.random.uniform(-2, 2, size=N)
        s.particles.image[:] = numpy.random.randint(-8, 8, size=(N, 3))
        s.particles.orientation[:] = numpy.random.uniform(-1, 1, size=(N, 4))
        s.particles.moment_inertia[:] = numpy.random.uniform(1, 5, size=(N, 3))
        s.particles.angmom[:] = numpy.random.uniform(-1, 1, size=(N, 4))
        s.particles.types = ['A', 'B', 'C', 'D']

        s.bonds.N = N - 1
        for i in range(s.bonds.N):
            s.bonds.group[i, :] = [i, i + 1]

        s.bonds.typeid[:] = numpy.random.randint(0, 3, size=s.bonds.N)
        s.bonds.types = ['bondA', 'bondB', 'bondC', 'bondD']

        s.angles.N = N - 2
        for i in range(s.angles.N):
            s.angles.group[i, :] = [i, i + 1, i + 2]

        s.angles.typeid[:] = numpy.random.randint(0, 3, size=s.angles.N)
        s.angles.types = ['angleA', 'angleB', 'angleC', 'angleD']

        s.dihedrals.N = N - 3
        for i in range(s.dihedrals.N):
            s.dihedrals.group[i, :] = [i, i + 1, i + 2, i + 3]

        s.dihedrals.typeid[:] = numpy.random.randint(0, 3, size=s.dihedrals.N)
        s.dihedrals.types = ['dihedralA', 'dihedralB', 'dihedralC', 'dihedralD']

        s.impropers.N = N - 3
        for i in range(s.impropers.N):
            s.impropers.group[i, :] = [i, i + 1, i + 2, i + 3]

        s.impropers.typeid[:] = numpy.random.randint(0, 3, size=s.impropers.N)
        s.impropers.types = ['improperA', 'improperB', 'improperC', 'improperD']

        s.pairs.N = N - 1
        for i in range(s.pairs.N):
            s.pairs.group[i, :] = [i, i + 1]

        s.pairs.typeid[:] = numpy.random.randint(0, 3, size=s.pairs.N)
        s.pairs.types = ['pairA', 'pairB', 'pairC', 'pairD']

        s.constraints.N = N - 1
        for i in range(s.constraints.N):
            s.constraints.group[i, :] = [i, i + 1]

        s.constraints.value[:] = numpy.random.uniform(1,
                                                      10,
                                                      size=s.constraints.N)

    return s

def assert_snapshots_equal(s1, s2):
    if s1.communicator.rank == 0:
        numpy.testing.assert_allclose(s1.configuration.box,
                                      s2.configuration.box)
        numpy.testing.assert_allclose(s1.configuration.dimensions,
                                      s2.configuration.dimensions)

        assert s1.particles.N == s2.particles.N
        assert s1.particles.types == s2.particles.types
        numpy.testing.assert_allclose(s1.particles.position,
                                      s2.particles.position)
        numpy.testing.assert_allclose(s1.particles.velocity,
                                      s2.particles.velocity)
        numpy.testing.assert_allclose(s1.particles.acceleration,
                                      s2.particles.acceleration)
        numpy.testing.assert_equal(s1.particles.typeid, s2.particles.typeid)
        numpy.testing.assert_allclose(s1.particles.mass, s2.particles.mass)
        numpy.testing.assert_allclose(s1.particles.charge, s2.particles.charge)
        numpy.testing.assert_allclose(s1.particles.diameter,
                                      s2.particles.diameter)
        numpy.testing.assert_equal(s1.particles.image, s2.particles.image)
        numpy.testing.assert_equal(s1.particles.body, s2.particles.body)
        numpy.testing.assert_allclose(s1.particles.orientation,
                                      s2.particles.orientation)
        numpy.testing.assert_allclose(s1.particles.moment_inertia,
                                      s2.particles.moment_inertia)
        numpy.testing.assert_allclose(s1.particles.angmom, s2.particles.angmom)
        numpy.testing.assert_allclose(s1.particles.diameter,
                                      s2.particles.diameter)

        assert s1.bonds.N == s2.bonds.N
        assert s1.bonds.types == s2.bonds.types
        numpy.testing.assert_equal(s1.bonds.typeid, s2.bonds.typeid)
        numpy.testing.assert_equal(s1.bonds.group, s2.bonds.group)

        assert s1.angles.N == s2.angles.N
        assert s1.angles.types == s2.angles.types
        numpy.testing.assert_equal(s1.angles.typeid, s2.angles.typeid)
        numpy.testing.assert_equal(s1.angles.group, s2.angles.group)

        assert s1.dihedrals.N == s2.dihedrals.N
        assert s1.dihedrals.types == s2.dihedrals.types
        numpy.testing.assert_equal(s1.dihedrals.typeid, s2.dihedrals.typeid)
        numpy.testing.assert_equal(s1.dihedrals.group, s2.dihedrals.group)

        assert s1.impropers.N == s2.impropers.N
        assert s1.impropers.types == s2.impropers.types
        numpy.testing.assert_equal(s1.impropers.typeid, s2.impropers.typeid)
        numpy.testing.assert_equal(s1.impropers.group, s2.impropers.group)

        assert s1.pairs.N == s2.pairs.N
        assert s1.pairs.types == s2.pairs.types
        numpy.testing.assert_equal(s1.pairs.typeid, s2.pairs.typeid)
        numpy.testing.assert_equal(s1.pairs.group, s2.pairs.group)

        assert s1.constraints.N == s2.constraints.N
        numpy.testing.assert_allclose(s1.constraints.value,
                                      s2.constraints.value)
        numpy.testing.assert_equal(s1.constraints.group, s2.constraints.group)

def test_create_from_snapshot(simulation_factory, snap):
    sim = simulation_factory(snap)

    if snap.communicator.rank == 0:
        assert sim.state.particle_types == snap.particles.types
        assert sim.state.bond_types == snap.bonds.types
        assert sim.state.angle_types == snap.angles.types
        assert sim.state.dihedral_types == snap.dihedrals.types
        assert sim.state.improper_types == snap.impropers.types
        assert sim.state.special_pair_types == snap.pairs.types
        # TODO: test box, dimensions
        # assert
        # TODO: what other state properties should be accessible and valid here?

def test_get_snapshot(simulation_factory, snap):
    sim = simulation_factory()
    sim.create_state_from_snapshot(snap)

    snap2 = sim.state.snapshot
    assert_snapshots_equal(snap, snap2)

def test_modify_snapshot(simulation_factory, snap):
    sim = simulation_factory()
    sim.create_state_from_snapshot(snap)

    if snap.communicator.rank == 0:
        snap.particles.N = snap.particles.N // 2
        snap.bonds.N = snap.bonds.N // 4
        snap.angles.N = snap.angles.N // 4
        snap.dihedrals.N = snap.dihedrals.N // 4
        snap.impropers.N = snap.impropers.N // 4
        snap.pairs.N = snap.pairs.N // 4
        snap.constraints.N = snap.constraints.N // 4

    sim.state.snapshot = snap

    snap2 = sim.state.snapshot
    assert_snapshots_equal(snap, snap2)

def test_thermalize_particle_velocity(simulation_factory,
                                      lattice_snapshot_factory):
    snap = lattice_snapshot_factory()
    sim = simulation_factory(snap)
    sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=1.5)

    snapshot = sim.state.snapshot
    if snapshot.communicator.rank == 0:
        v = snapshot.particles.velocity[:]
        m = snapshot.particles.mass[:]
        p = m * v.T
        p_com = numpy.mean(p, axis=1)

        numpy.testing.assert_allclose(p_com, [0, 0, 0], atol=1e-14)

        K = numpy.sum(1 / 2 * m * (v[:, 0]**2 + v[:, 1]**2 + v[:, 2]**2))
        # check that K is somewhat close to the target - the fluctuations are
        # too large for an allclose check.
        expected_K = (3 * snap.particles.N - 3) / 2 * 1.5
        assert K > expected_K * 3 / 4 and K < expected_K * 4 / 3

def test_thermalize_angular_momentum(simulation_factory,
                                     lattice_snapshot_factory):
    snap = lattice_snapshot_factory()
    I = [1, 2, 3]  # noqa: E741 - allow ambiguous variable name

    if snap.communicator.rank == 0:
        snap.particles.moment_inertia[:] = I

    sim = simulation_factory(snap)
    sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=1.5)

    snapshot = sim.state.snapshot
    if snapshot.communicator.rank == 0:
        # Note: this conversion assumes that all particles have (1, 0, 0, 0)
        # orientations.
        L = snapshot.particles.angmom[:, 1:4] / 2

        K = numpy.sum(
            1 / 2 * (L[:, 0]**2 / I[0] + L[:, 1]**2 / I[1] + L[:, 2]**2 / I[2]))
        # check that K is somewhat close to the target - the fluctuations are
        # too large for an allclose check.
        expected_K = (3 * snap.particles.N) / 2 * 1.5
        assert K > expected_K * 3 / 4 and K < expected_K * 4 / 3

max_stars_repo_name: XT-Lee/hoomd-blue

max_issues_repo_name: XT-Lee/hoomd-blue

max_forks_repo_name: XT-Lee/hoomd-blue


text:
"""
Gaussian Mixture Models
"""

# Author: Ron Weiss <ronweiss@gmail.com>
#         Fabian Pedregosa <fabian.pedregosa@inria.fr>
#

import numpy as np

from sklearn.base import BaseEstimator
from sklearn import cluster

def logsum(A, axis=None):
    """Computes the sum of A assuming A is in the log domain.

    Returns log(sum(exp(A), axis)) while minimizing the possibility of
    over/underflow.
    """
    Amax = A.max(axis)
    if axis and A.ndim > 1:
        shape = list(A.shape)
        shape[axis] = 1
        Amax.shape = shape
    Asum = np.log(np.sum(np.exp(A - Amax), axis))
    Asum += Amax.reshape(Asum.shape)
    if axis:
        # Look out for underflow.
        Asum[np.isnan(Asum)] = - np.Inf
    return Asum

# TODO: this lacks a docstring
def normalize(A, axis=None):
    A += np.finfo(float).eps
    Asum = A.sum(axis)
    if axis and A.ndim > 1:
        # Make sure we don't divide by zero.
        Asum[Asum == 0] = 1
        shape = list(A.shape)
        shape[axis] = 1
        Asum.shape = shape
    return A / Asum

def lmvnpdf(obs, means, covars, cvtype='diag'):
    """Compute the log probability under a multivariate Gaussian distribution.

    Parameters
    ----------
    obs : array_like, shape (O, D)
        List of D-dimensional data points.  Each row corresponds to a
        single data point.
    means : array_like, shape (C, D)
        List of D-dimensional mean vectors for C Gaussians.  Each row
        corresponds to a single mean vector.
    covars : array_like
        List of C covariance parameters for each Gaussian.  The shape
        depends on `cvtype`:
            (C,)      if 'spherical',
            (D, D)    if 'tied',
            (C, D)    if 'diag',
            (C, D, D) if 'full'
    cvtype : string
        Type of the covariance parameters.  Must be one of
        'spherical', 'tied', 'diag', 'full'.  Defaults to 'diag'.

    Returns
    -------
    lpr : array_like, shape (O, C)
        Array containing the log probabilities of each data point in
        `observation` under each of the C multivariate Gaussian distributions.
    """
    lmvnpdf_dict = {'spherical': _lmvnpdfspherical,
                    'tied': _lmvnpdftied,
                    'diag': _lmvnpdfdiag,
                    'full': _lmvnpdffull}
    return lmvnpdf_dict[cvtype](obs, means, covars)

def sample_gaussian(mean, covar, cvtype='diag', n_samples=1):
    """Generate random samples from a Gaussian distribution.

    Parameters
    ----------
    mean : array_like, shape (n_features,)
        Mean of the distribution.

    covars : array_like, optional
        Covariance of the distribution.  The shape depends on `cvtype`:
            scalar  if 'spherical',
            (D)     if 'diag',
            (D, D)  if 'tied', or 'full'

    cvtype : string, optional
        Type of the covariance parameters.  Must be one of
        'spherical', 'tied', 'diag', 'full'.  Defaults to 'diag'.

    n_samples : int, optional
        Number of samples to generate. Defaults to 1.

    Returns
    -------
    observation : array, shape (n_features, n_samples)
        Randomly generated sample
    """
    n_dim = len(mean)
    rand = np.random.randn(n_dim, n_samples)
    if n_samples == 1:
        rand.shape = (n_dim,)

    if cvtype == 'spherical':
        rand *= np.sqrt(covar)
    elif cvtype == 'diag':
        rand = np.dot(np.diag(np.sqrt(covar)), rand)
    else:
        from scipy import linalg
        U, s, V = linalg.svd(covar)
        sqrtS = np.diag(np.sqrt(s))
        sqrt_covar = np.dot(U, np.dot(sqrtS, V))
        rand = np.dot(sqrt_covar, rand)

    return (rand.T + mean).T

class GMM(BaseEstimator):
    """Gaussian Mixture Model

    Representation of a Gaussian mixture model probability distribution.
    This class allows for easy evaluation of, sampling from, and
    maximum-likelihood estimation of the parameters of a GMM distribution.

    Initializes parameters such that every mixture component has zero
    mean and identity covariance.

    Parameters
    ----------
    n_states : int, optional
        Number of mixture components. Defaults to 1.

    cvtype : string (read-only), optional
        String describing the type of covariance parameters to
        use.  Must be one of 'spherical', 'tied', 'diag', 'full'.
        Defaults to 'diag'.

    Attributes
    ----------
    cvtype : string (read-only)
        String describing the type of covariance parameters used by
        the GMM.  Must be one of 'spherical', 'tied', 'diag', 'full'.
    n_features : int
        Dimensionality of the Gaussians.
    n_states : int (read-only)
        Number of mixture components.
    weights : array, shape (`n_states`,)
        Mixing weights for each mixture component.
    means : array, shape (`n_states`, `n_features`)
        Mean parameters for each mixture component.
    covars : array
        Covariance parameters for each mixture component.  The shape
        depends on `cvtype`:
            (`n_states`,)                             if 'spherical',
            (`n_features`, `n_features`)              if 'tied',
            (`n_states`, `n_features`)                if 'diag',
            (`n_states`, `n_features`, `n_features`)  if 'full'
    converged_ : bool
        True when convergence was reached in fit(), False
        otherwise.

    Methods
    -------
    decode(X)
        Find most likely mixture components for each point in `X`.
    eval(X)
        Compute the log likelihood of `X` under the model and the
        posterior distribution over mixture components.
    fit(X)
        Estimate model parameters from `X` using the EM algorithm.
    predict(X)
        Like decode, find most likely mixtures components for each
        observation in `X`.
    rvs(n=1)
        Generate `n` samples from the model.
    score(X)
        Compute the log likelihood of `X` under the model.

    Examples
    --------
    >>> import numpy as numpy
    >>> g = mixture.GMM(n_states=2)

    >>> # Generate random observations with two modes centered on 0
    >>> # and 10 to use for training.
    >>> numpy.random.seed(0)
    >>> observation = numpy.concatenate((numpy.random.randn(100, 1),
    ...                       10 + numpy.random.randn(300, 1)))
    >>> g.fit(observation)
    GMM(cvtype='diag', n_states=2)
    >>> g.weights
    array([ 0.25,  0.75])
    >>> g.means
    array([[ 0.05980802],
           [ 9.94199467]])
    >>> g.covars
    [array([[ 1.01682662]]), array([[ 0.96080513]])]
    >>> numpy.round(g.weights, 2)
    array([ 0.25,  0.75])
    >>> numpy.round(g.means, 2)
    array([[ 0.06],
           [ 9.94]])
    >>> numpy.round(g.covars, 2)
    ... #doctest: +NORMALIZE_WHITESPACE
    array([[[ 1.02]],
           [[ 0.96]]])
    >>> g.predict([[0], [2], [9], [10]])
    array([0, 0, 1, 1])
    >>> numpy.round(g.score([[0], [2], [9], [10]]), 2)
    array([-2.32, -4.16, -1.65, -1.19])

    >>> # Refit the model on new data (initial parameters remain the
    >>> # same), this time with an even split between the two modes.
    >>> g.fit(20 * [[0]] +  20 * [[10]])
    GMM(cvtype='diag', n_states=2)
    >>> numpy.round(g.weights, 2)
    array([ 0.5,  0.5])
    """

    def __init__(self, n_states=1, cvtype='diag'):
        self._n_states = n_states
        self._cvtype = cvtype

        if not cvtype in ['spherical', 'tied', 'diag', 'full']:
            raise ValueError('bad cvtype')

        self.weights = np.ones(self._n_states) / self._n_states

        # flag to indicate exit status of fit() method: converged (True) or
        # n_iter reached (False)
        self.converged_ = False

    # Read-only properties.
    @property
    def cvtype(self):
        """Covariance type of the model.

        Must be one of 'spherical', 'tied', 'diag', 'full'.
        """
        return self._cvtype

    @property
    def n_states(self):
        """Number of mixture components in the model."""
        return self._n_states

    def _get_covars(self):
        """Return covars as a full matrix."""
        if self.cvtype == 'full':
            return self._covars
        elif self.cvtype == 'diag':
            return [np.diag(cov) for cov in self._covars]
        elif self.cvtype == 'tied':
            return [self._covars] * self._n_states
        elif self.cvtype == 'spherical':
            return [np.eye(self.n_features) * f for f in self._covars]

    def _set_covars(self, covars):
        covars = np.asanyarray(covars)
        _validate_covars(covars, self._cvtype, self._n_states, self.n_features)
        self._covars = covars

    covars = property(_get_covars, _set_covars)

    def _get_means(self):
        """Mean parameters for each mixture component."""
        return self._means

    def _set_means(self, means):
        means = np.asarray(means)
        if hasattr(self, 'n_features') and \
               means.shape != (self._n_states, self.n_features):
            raise ValueError('means must have shape (n_states, n_features)')
        self._means = means.copy()
        self.n_features = self._means.shape[1]

    means = property(_get_means, _set_means)

    def _get_weights(self):
        """Mixing weights for each mixture component."""
        return np.exp(self._log_weights)

    def _set_weights(self, weights):
        if len(weights) != self._n_states:
            raise ValueError('weights must have length n_states')
        if not np.allclose(np.sum(weights), 1.0):
            raise ValueError('weights must sum to 1.0')

        self._log_weights = np.log(np.asarray(weights).copy())

    weights = property(_get_weights, _set_weights)

    def eval(self, obs):
        """Evaluate the model on data

        Compute the log probability of `observation` under the model and
        return the posterior distribution (responsibilities) of each
        mixture component for each element of `observation`.

        Parameters
        ----------
        obs : array_like, shape (n_samples, n_features)
            List of n_features-dimensional data points.  Each row
            corresponds to a single data point.

        Returns
        -------
        logprob : array_like, shape (n_samples,)
            Log probabilities of each data point in `observation`
        posteriors: array_like, shape (n_samples, n_states)
            Posterior probabilities of each mixture component for each
            observation
        """
        obs = np.asanyarray(obs)
        lpr = (lmvnpdf(obs, self._means, self._covars, self._cvtype)
               + self._log_weights)
        logprob = logsum(lpr, axis=1)
        posteriors = np.exp(lpr - logprob[:, np.newaxis])
        return logprob, posteriors

    def score(self, obs):
        """Compute the log probability under the model.

        Parameters
        ----------
        obs : array_like, shape (n_samples, n_features)
            List of n_features-dimensional data points.  Each row
            corresponds to a single data point.

        Returns
        -------
        logprob : array_like, shape (n_samples,)
            Log probabilities of each data point in `observation`
        """
        logprob, posteriors = self.eval(obs)
        return logprob

    def decode(self, obs):
        """Find most likely mixture components for each point in `observation`.

        Parameters
        ----------
        obs : array_like, shape (n, n_features)
            List of n_features-dimensional data points.  Each row
            corresponds to a single data point.

        Returns
        -------
        logprobs : array_like, shape (n_samples,)
            Log probability of each point in `observation` under the model.
        components : array_like, shape (n_samples,)
            Index of the most likelihod mixture components for each observation
        """
        logprob, posteriors = self.eval(obs)
        return logprob, posteriors.argmax(axis=1)

    def predict(self, X):
        """Predict label for data.

        Parameters
        ----------
        X : array-like, shape = [n_samples, n_features]

        Returns
        -------
        C : array, shape = (n_samples,)
        """
        logprob, components = self.decode(X)
        return components

    def predict_proba(self, X):
        """Predict posterior probability of data under each Gaussian
        in the model.

        Parameters
        ----------
        X : array-like, shape = [n_samples, n_features]

        Returns
        -------
        T : array-like, shape = (n_samples, n_states)
            Returns the probability of the sample for each Gaussian
            (state) in the model.
        """
        logprob, posteriors = self.eval(X)
        return posteriors

    def rvs(self, n_samples=1):
        """Generate random samples from the model.

        Parameters
        ----------
        n_samples : int, optional
            Number of samples to generate. Defaults to 1.

        Returns
        -------
        observation : array_like, shape (n_samples, n_features)
            List of samples
        """
        weight_pdf = self.weights
        weight_cdf = np.cumsum(weight_pdf)

        obs = np.empty((n_samples, self.n_features))
        rand = np.random.rand(n_samples)
        # decide which component to use for each sample
        comps = weight_cdf.searchsorted(rand)
        # for each component, generate all needed samples
        for comp in range(self._n_states):
            # occurrences of current component in observation
            comp_in_obs = (comp == comps)
            # number of those occurrences
            num_comp_in_obs = comp_in_obs.sum()
            if num_comp_in_obs > 0:
                if self._cvtype == 'tied':
                    cv = self._covars
                else:
                    cv = self._covars[comp]
                obs[comp_in_obs] = sample_gaussian(
                    self._means[comp], cv, self._cvtype, num_comp_in_obs).T
        return obs

    def fit(self, X, n_iter=10, min_covar=1e-3, thresh=1e-2, params='wmc',
            init_params='wmc'):
        """Estimate model parameters with the expectation-maximization
        algorithm.

        A initialization step is performed before entering the em
        algorithm. If you want to avoid this step, set the keyword
        argument init_params to the empty string ''. Likewise, if you
        would like just to do an initialization, call this method with
        n_iter=0.

        Parameters
        ----------
        X : array_like, shape (n, n_features)
            List of n_features-dimensional data points.  Each row
            corresponds to a single data point.
        n_iter : int, optional
            Number of EM iterations to perform.
        min_covar : float, optional
            Floor on the diagonal of the covariance matrix to prevent
            overfitting.  Defaults to 1e-3.
        thresh : float, optional
            Convergence threshold.
        params : string, optional
            Controls which parameters are updated in the training
            process.  Can contain any combination of 'w' for weights,
            'm' for means, and 'c' for covars.  Defaults to 'wmc'.
        init_params : string, optional
            Controls which parameters are updated in the initialization
            process.  Can contain any combination of 'w' for weights,
            'm' for means, and 'c' for covars.  Defaults to 'wmc'.
        """

        ## initialization step

        X = np.asanyarray(X)

        if hasattr(self, 'n_features') and self.n_features != X.shape[1]:
            raise ValueError('Unexpected number of dimensions, got %s but '
                             'expected %s' % (X.shape[1], self.n_features))

        self.n_features = X.shape[1]

        if 'm' in init_params:
            self._means = cluster.KMeans(
                k=self._n_states).fit(X).cluster_centers_
        elif not hasattr(self, 'means'):
                self._means = np.zeros((self.n_states, self.n_features))

        if 'w' in init_params or not hasattr(self, 'weights'):
            self.weights = np.tile(1.0 / self._n_states, self._n_states)

        if 'c' in init_params:
            cv = np.cov(X.T)
            if not cv.shape:
                cv.shape = (1, 1)
            self._covars = _distribute_covar_matrix_to_match_cvtype(
                cv, self._cvtype, self._n_states)
        elif not hasattr(self, 'covars'):
                self.covars = _distribute_covar_matrix_to_match_cvtype(
                    np.eye(self.n_features), self.cvtype, self.n_states)

        # EM algorithm
        logprob = []
        # reset self.converged_ to False
        self.converged_ = False
        for i in range(n_iter):
            # Expectation step
            curr_logprob, posteriors = self.eval(X)
            logprob.append(curr_logprob.sum())

            # Check for convergence.
            if i > 0 and abs(logprob[-1] - logprob[-2]) < thresh:
                self.converged_ = True
                break

            # Maximization step
            self._do_mstep(X, posteriors, params, min_covar)

        return self

    def _do_mstep(self, X, posteriors, params, min_covar=0):
        w = posteriors.sum(axis=0)
        avg_obs = np.dot(posteriors.T, X)
        norm = 1.0 / (w[:, np.newaxis] + 1e-200)

        if 'w' in params:
            self._log_weights = np.log(w / w.sum())
        if 'm' in params:
            self._means = avg_obs * norm
        if 'c' in params:
            covar_mstep_func = _covar_mstep_funcs[self._cvtype]
            self._covars = covar_mstep_func(self, X, posteriors,
                                            avg_obs, norm, min_covar)

        return w

##
## some helper routines
##

def _lmvnpdfdiag(obs, means=0.0, covars=1.0):
    n_obs, n_dim = obs.shape
    # (x-y).T A (x-y) = x.T A x - 2x.T A y + y.T A y
    #lpr = -0.5 * (numpy.tile((numpy.sum((means**2) / covars, 1)
    #                  + numpy.sum(numpy.log(covars), 1))[numpy.newaxis,:], (n_obs,1))
    lpr = -0.5 * (n_dim * np.log(2 * np.pi) + np.sum(np.log(covars), 1)
                  + np.sum((means ** 2) / covars, 1)
                  - 2 * np.dot(obs, (means / covars).T)
                  + np.dot(obs ** 2, (1.0 / covars).T))
    return lpr

def _lmvnpdfspherical(obs, means=0.0, covars=1.0):
    cv = covars.copy()
    if covars.ndim == 1:
        cv = cv[:, np.newaxis]
    return _lmvnpdfdiag(obs, means, np.tile(cv, (1, obs.shape[-1])))

def _lmvnpdftied(obs, means, covars):
    from scipy import linalg
    n_obs, n_dim = obs.shape
    # (x-y).T A (x-y) = x.T A x - 2x.T A y + y.T A y
    icv = linalg.pinv(covars)
    lpr = -0.5 * (n_dim * np.log(2 * np.pi) + np.log(linalg.det(covars))
                  + np.sum(obs * np.dot(obs, icv), 1)[:, np.newaxis]
                  - 2 * np.dot(np.dot(obs, icv), means.T)
                  + np.sum(means * np.dot(means, icv), 1))
    return lpr

def _lmvnpdffull(obs, means, covars):
    """
    Log probability for full covariance matrices.
    """
    from scipy import linalg
    import itertools
    if hasattr(linalg, 'solve_triangular'):
        # only in scipy since 0.9
        solve_triangular = linalg.solve_triangular
    else:
        # slower, but works
        solve_triangular = linalg.solve
    n_obs, n_dim = obs.shape
    nmix = len(means)
    log_prob = np.empty((n_obs, nmix))
    for c, (mu, cv) in enumerate(itertools.izip(means, covars)):
        cv_chol = linalg.cholesky(cv, lower=True)
        cv_log_det = 2 * np.sum(np.log(np.diagonal(cv_chol)))
        cv_sol = solve_triangular(cv_chol, (obs - mu).T, lower=True).T
        log_prob[:, c] = -.5 * (np.sum(cv_sol ** 2, axis=1) + \
                           n_dim * np.log(2 * np.pi) + cv_log_det)

    return log_prob

def _validate_covars(covars, cvtype, nmix, n_dim):
    from scipy import linalg
    if cvtype == 'spherical':
        if len(covars) != nmix:
            raise ValueError("'spherical' covars must have length nmix")
        elif np.any(covars <= 0):
            raise ValueError("'spherical' covars must be non-negative")
    elif cvtype == 'tied':
        if covars.shape != (n_dim, n_dim):
            raise ValueError("'tied' covars must have shape (n_dim, n_dim)")
        elif (not np.allclose(covars, covars.T)
              or np.any(linalg.eigvalsh(covars) <= 0)):
            raise ValueError("'tied' covars must be symmetric, "
                             "positive-definite")
    elif cvtype == 'diag':
        if covars.shape != (nmix, n_dim):
            raise ValueError("'diag' covars must have shape (nmix, n_dim)")
        elif np.any(covars <= 0):
            raise ValueError("'diag' covars must be non-negative")
    elif cvtype == 'full':
        if covars.shape != (nmix, n_dim, n_dim):
            raise ValueError("'full' covars must have shape "
                             "(nmix, n_dim, n_dim)")
        for n, cv in enumerate(covars):
            if (not np.allclose(cv, cv.T)
                or np.any(linalg.eigvalsh(cv) <= 0)):
                raise ValueError("component %d of 'full' covars must be "
                                 "symmetric, positive-definite" % n)

def _distribute_covar_matrix_to_match_cvtype(tiedcv, cvtype, n_states):
    if cvtype == 'spherical':
        cv = np.tile(np.diag(tiedcv).mean(), n_states)
    elif cvtype == 'tied':
        cv = tiedcv
    elif cvtype == 'diag':
        cv = np.tile(np.diag(tiedcv), (n_states, 1))
    elif cvtype == 'full':
        cv = np.tile(tiedcv, (n_states, 1, 1))
    else:
        raise (ValueError,
               "cvtype must be one of 'spherical', 'tied', 'diag', 'full'")
    return cv

def _covar_mstep_diag(gmm, obs, posteriors, avg_obs, norm, min_covar):
    # For column vectors:
    # covars_c = average((observation(t) - means_c) (observation(t) - means_c).T,
    #                    weights_c)
    # (observation(t) - means_c) (observation(t) - means_c).T
    #     = observation(t) observation(t).T - 2 observation(t) means_c.T + means_c means_c.T
    #
    # But everything here is a row vector, so all of the
    # above needs to be transposed.
    avg_obs2 = np.dot(posteriors.T, obs * obs) * norm
    avg_means2 = gmm._means ** 2
    avg_obs_means = gmm._means * avg_obs * norm
    return avg_obs2 - 2 * avg_obs_means + avg_means2 + min_covar

def _covar_mstep_spherical(*args):
    return _covar_mstep_diag(*args).mean(axis=1)

def _covar_mstep_full(gmm, obs, posteriors, avg_obs, norm, min_covar):
    # Eq. 12 from K. Murphy, "Fitting a Conditional Linear Gaussian
    # Distribution"
    cv = np.empty((gmm._n_states, gmm.n_features, gmm.n_features))
    for c in range(gmm._n_states):
        post = posteriors[:, c]
        avg_cv = np.dot(post * obs.T, obs) / post.sum()
        mu = gmm._means[c][np.newaxis]
        cv[c] = (avg_cv - np.dot(mu.T, mu)
                 + min_covar * np.eye(gmm.n_features))
    return cv

def _covar_mstep_tied2(*args):
    return _covar_mstep_full(*args).mean(axis=0)

def _covar_mstep_tied(gmm, obs, posteriors, avg_obs, norm, min_covar):
    #print "THIS IS BROKEN"
    # Eq. 15 from K. Murphy, "Fitting a Conditional Linear Gaussian
    avg_obs2 = np.dot(obs.T, obs)
    avg_means2 = np.dot(gmm._means.T, gmm._means)
    return (avg_obs2 - avg_means2 + min_covar * np.eye(gmm.n_features))

def _covar_mstep_slow(gmm, obs, posteriors, avg_obs, norm, min_covar):
    w = posteriors.sum(axis=0)
    covars = np.zeros(gmm._covars.shape)
    for c in range(gmm._n_states):
        mu = gmm._means[c]
        #cv = numpy.dot(mu.T, mu)
        avg_obs2 = np.zeros((gmm.n_features, gmm.n_features))
        for t, o in enumerate(obs):
            avg_obs2 += posteriors[t, c] * np.outer(o, o)
        cv = (avg_obs2 / w[c]
              - 2 * np.outer(avg_obs[c] / w[c], mu)
              + np.outer(mu, mu)
              + min_covar * np.eye(gmm.n_features))
        if gmm.cvtype == 'spherical':
            covars[c] = np.diag(cv).mean()
        elif gmm.cvtype == 'diag':
            covars[c] = np.diag(cv)
        elif gmm.cvtype == 'full':
            covars[c] = cv
        elif gmm.cvtype == 'tied':
            covars += cv / gmm._n_states
    return covars

_covar_mstep_funcs = {'spherical': _covar_mstep_spherical,
                      'diag': _covar_mstep_diag,
                      #'tied': _covar_mstep_tied,
                      'full': _covar_mstep_full,
                      'tied': _covar_mstep_slow,
                      }


text:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import sys
sys.path.insert(0, './HydroBayesOpt')
from bayes_regression import BayesianRegression
from gp import GaussianProcesses
from matplotlib.colors import rgb2hex
np.random.seed(3)

def run_bayesian_regression_experiment():
	# Two models
	function_linear = lambda x: np.array([[1, xx] for xx in x])
	function_polynomial = lambda x: np.array([[xx**4, 4.*xx**3, 8.*xx, 1] for xx in x])

	n_features = 4
	function = function_polynomial

	# Create priors
	prior_V = np.eye(n_features, dtype=float)
	prior_mu = np.zeros(n_features)

	# Create dataset
	N = [3, 5, 15, 50]
	bounds = [-4., 2.]
	var_error = 5.

	# Generate random data set
	x0 = np.random.uniform(bounds[0], bounds[1], N[-1])
	y0 = np.sum(function(x0), axis=1)
	y0 += np.random.normal(0, var_error, N[-1])

	# Create regression object
	regression = BayesianRegression(var_error, prior_mu, prior_V, function=function_polynomial)

	# Add data to regression model in stages and sample model from parameter posterior
	fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)
	n_previous = 0
	for ax, n in zip(axes.ravel(), N):
		regression.add_observations(x0[n_previous:n], y0[n_previous:n])
		n_previous = n
		regression.plot_sample_models_1d(ax, 5)
		ax.set_title('N={}'.format(n))

	# Create new regression object
	regression = BayesianRegression(var_error, prior_mu, prior_V, function=function_polynomial)

	# Add data to regression model in stages and plot model mean and error variance
	plt.figure(1)
	fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)
	n_previous = 0
	for ax, n in zip(axes.ravel(), N):
		regression.add_observations(x0[n_previous:n], y0[n_previous:n])
		n_previous = n
		regression.plot_predictive_model(ax)
		ax.set_title('N={}'.format(n))

	# plot parameter distirbution after all data is added
	plt.figure(2)
	regression.plot_parameter_distributions()

	plt.show()

def run_gaussian_processes_experiment():
	boundaries = [-5., 5.01]
	gp = GaussianProcesses(0., boundaries, kernel='rbf', kernel_params=1.)
	n_samples = 5
	dx = 2

	X = np.array(np.arange(boundaries[0], boundaries[1], np.ptp(boundaries) / n_samples)) + 0.5 + np.random.uniform(-0.2, 0.2, n_samples)
	y = np.array(np.random.uniform(-1, 1, n_samples))

	gp.add_data(X, y)
	samples, x_samples = gp.sample_functions(3, (boundaries[1] - boundaries[0]) / dx)

	plt.figure(1)

	# plt.figure(2)
	cmap = cm.get_cmap('Set2')
	x_star = np.arange(boundaries[0], boundaries[1], dx)
	y_star, sigma_star = gp.predict(x_star)
	# sigma_star = 2. * sigma_star ** .5
	plt.fill_between(x_star, y_star + sigma_star, y_star - sigma_star, facecolor='lightgrey', label='$\mu \pm \sigma$')
	plt.plot(x_star, y_star, c='r', label='$\mu$')
	plt.plot(x_star, y_star + sigma_star, c='darkgrey', alpha=0.4)
	plt.plot(x_star, y_star - sigma_star, c='darkgrey', alpha=0.4)
	plt.scatter(X, y, c=rgb2hex(cmap(2)), s=20)
	plt.ylim([-1.5, 1.5])
	plt.xlim(boundaries)
	plt.legend()
	plt.xlabel('x')

	plt.figure(2)
	plt.fill_between(x_star, y_star + sigma_star, y_star - sigma_star, facecolor='lightgrey', label='$\mu \pm \sigma$')
	plt.plot(x_star, y_star, c='darkgrey')
	plt.plot(x_star, y_star + sigma_star, c='darkgrey', alpha=0.4)
	plt.plot(x_star, y_star - sigma_star, c='darkgrey', alpha=0.4)
	colors = [cmap(i) for i in range(3)]
	for s, c in zip(samples, colors):
		plt.plot(x_samples, s, c=c)#, alpha=0.2)
	plt.ylim([-1.5, 1.5])
	plt.xlim(boundaries)
	plt.legend()
	plt.xlabel('y')

	plt.show()

# run_bayesian_regression_experiment()
run_gaussian_processes_experiment()

max_stars_repo_name: bernardoct/BayesianRegressionAndGPs

max_issues_repo_name: bernardoct/BayesianRegressionAndGPs

max_forks_repo_name: bernardoct/BayesianRegressionAndGPs


text:
"""
Very simple training script for a mask estimator.
Saves checkpoints and summaries to $STORAGE_ROOT/speech_enhancement/simple_mask_estimator_{id}
may be called with:
python -m padertorch.contrib.examples.speech_enhancement.simple_mask_estimator.train with database_json=/path/to/json
"""

from pathlib import Path

import os
import numpy as np
import paderbox as pb
import padertorch as pt
from lazy_dataset.database import JsonDatabase
from pb_bss.extraction.mask_module import biased_binary_mask
from sacred import Experiment, observers

from .model import SimpleMaskEstimator

ex = Experiment('Train Simple Mask Estimator')

@ex.config
def config():
    storage_dir = None
    if storage_dir is None:
        storage_dir = pt.io.get_new_storage_dir(
            'speech_enhancement', prefix='simple_mask_estimator')
    database_json = None
    if database_json is None:
        if 'NT_DATABASE_JSONS_DIR' in os.environ:
            database_json = Path(
                os.environ['NT_DATABASE_JSONS_DIR']) / 'chime.json'
    assert database_json is not None, (
        'You have to specify a path to a json describing your database,'
        'use "with database_json=/Path/To/Json" as suffix to your call'
    )
    assert Path(database_json).exists(), database_json
    ex.observers.append(observers.FileStorageObserver(
        Path(storage_dir).expanduser().resolve() / 'sacred')
    )

def prepare_data(example):
    stft = pb.transform.STFT(shift=256, size=1024)
    net_input = dict()
    audio_data = dict()
    for key in ['observation', 'speech_image', 'noise_image']:
        audio_data[key] = stft(np.array([
            pb.io.load_audio(audio) for audio in example['audio_path'][key]]))
    net_input['observation_abs'] = np.abs(
        audio_data['observation']).astype(np.float32)
    target_mask, noise_mask = biased_binary_mask(np.stack(
        [audio_data['speech_image'], audio_data['noise_image']], axis=0
    ))
    net_input['speech_mask_target'] = target_mask.astype(np.float32)
    net_input['noise_mask_target'] = noise_mask.astype(np.float32)
    return net_input

def get_train_dataset(database: JsonDatabase):
    train_ds = database.get_dataset('tr05_simu')
    return (train_ds
            .map(prepare_data)
            .prefetch(num_workers=4, buffer_size=4))

def get_validation_dataset(database: JsonDatabase):
    # AudioReader is a specialized function to read audio organized
    # in a json as described in pb.database.database
    val_iterator = database.get_dataset('dt05_simu')
    return val_iterator.map(prepare_data) \
        .prefetch(num_workers=4, buffer_size=4)

@ex.command
def test_run(storage_dir, database_json):
    model = SimpleMaskEstimator(513)
    print(f'Simple training for the following model: {model}')
    database = JsonDatabase(database_json)
    train_dataset = get_train_dataset(database)
    validation_dataset = get_validation_dataset(database)
    trainer = pt.train.trainer.Trainer(
        model, storage_dir, optimizer=pt.train.optimizer.Adam(),
        stop_trigger=(int(1e5), 'iteration')
    )
    trainer.test_run(train_dataset, validation_dataset)

@ex.automain
def train(storage_dir, database_json):
    model = SimpleMaskEstimator(513)
    print(f'Simple training for the following model: {model}')
    database = JsonDatabase(database_json)
    train_dataset = get_train_dataset(database)
    validation_dataset = get_validation_dataset(database)
    trainer = pt.Trainer(model, storage_dir,
                         optimizer=pt.train.optimizer.Adam(),
                         stop_trigger=(int(1e5), 'iteration'))
    trainer.test_run(train_dataset, validation_dataset)
    trainer.register_validation_hook(
        validation_dataset, n_back_off=5, lr_update_factor=1 / 10,
        back_off_patience=1, early_stopping_patience=None)
    trainer.train(train_dataset)

max_stars_repo_name: sibange/padertorch

max_issues_repo_name: sibange/padertorch

max_forks_repo_name: sibange/padertorch