text:
import os
import sys
import scipy.io
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_moons, make_circles, make_classification
from sklearn.neural_network import MLPClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.svm import SVC
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.gaussian_process.kernels import RBF
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
from sklearn.metrics import *

from sklearn.model_selection import cross_val_predict
from sklearn.model_selection import cross_validate
from sklearn.model_selection import ShuffleSplit
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from sklearn.metrics import *
from sklearn import metrics

#from lib.Shared_Functions import *
import glob

#%          ################ FUNCTION   ###########################
def Get_model_performnace(y,y_predicted):

    if len(set(y ))==2:

        C = confusion_matrix(y,y_predicted)
    #    print('Confusion Matrix : \n', C)

        total1=sum(sum(C))
        #####from confusion matrix calculate accuracy
        accuracy=(C[0,0]+C[1,1])/total1
    #    print ('Accuracy : ', accuracy)

        sensitivity = C[0,0]/(C[0,0]+C[0,1])
    #    print('Sensitivity : ', sensitivity )

        specificity = C[1,1]/(C[1,0]+C[1,1])
    #    print('Specificity : ', specificity)

    #    # accuracy: (tp + tn) / (p + n)
    #    accuracy = accuracy_score(y, y_predicted)
    #    print('Accuracy: %f' % accuracy)
    #
    #
        # precision tp / (tp + fp)
        precision = precision_score(y, y_predicted)
    #    print('Precision: %f' % precision)
        # recall: tp / (tp + fn)
        recall = recall_score(y, y_predicted)
    #    print('Recall: %f' % recall)

        # f1: 2 tp / (2 tp + fp + fn)
        f1 = f1_score(y, y_predicted)
    #    print('F1 score: %f' % f1)
        fpr, tpr, thresholds = metrics.roc_curve(y, y_predicted)
        AUC=metrics.auc(fpr, tpr)
    #    print('AUC score: %f' % AUC)

    else:

        MC_performance=report2dict(classification_report(y, y_predicted))

        accuracy= round((y == y_predicted).mean().item(),2)
        specificity=-1
        recall=float(MC_performance['recall'][-1:])
        sensitivity=recall
        f1=float(MC_performance['f1-score'][-1:])
        precision=float(MC_performance['precision'][-1:])
        AUC=-1

#    print ('Accuracy : ', accuracy)

    return accuracy, sensitivity, specificity, precision, recall, f1, AUC #=Get_model_performnace(y_test,y_predicted)

#%%#################################################################################################
# input parameters
save_excel=1;
project_name='QuPWM'

#The classifier pipeline
names = ["Nearest Neighbors", "Linear SVM",
        "RBF SVM", "Decision Tree", "Random Forest", "Neural Net",# "AdaBoost",
        "Naive Bayes", "Quadratic Discriminant Analysis"]

classifiers = [KNeighborsClassifier(3), SVC(kernel="linear", C=0.025),
    SVC(gamma=2, C=1),
    DecisionTreeClassifier(max_depth=5),
    RandomForestClassifier(max_depth=5, n_estimators=10, max_features=1),
    MLPClassifier(alpha=1),#AdaBoostClassifier(),
    GaussianNB(),
    QuadraticDiscriminantAnalysis()]

#%% #################################################################################################
data_path='./mat/EA001EA002EA003EA004EA005EA006EA007EA008EA009/'
#data_path='./mat/'
#files = [f for f in glob.glob(data_path+"**/90/*.mat", recursive=True)]
#names, classifiers=["Nearest Neighbors"],[ KNeighborsClassifier(3)]
score=[];Clf_score=[];clf_k=[]
names_clf_L=[]
list_L=[80, 90,95, 100, 110]
for name, clf in zip(names, classifiers):
    for L in list_L:
        files = [f for f in glob.glob(data_path+"**/"+str(L) + "/*.mat", recursive=True)]
    #  ---------------------------- Train/Test Split  ----------------------------
        score_clf=[]; fold=0
        print('Cross Validation using:',name)

        for mat_filename in files:
            #print('File L=',L,' :',mat_filename)
            fold=fold+1

            #% The classifier pipelinLoad mat files
            loaded_mat_file = scipy.io.loadmat(mat_filename)
            X_train = loaded_mat_file['fPWM_train']
            X_test = loaded_mat_file['fPWM_test']
            y_train = loaded_mat_file['y_train'].ravel()
            y_test = loaded_mat_file['y_test'].ravel()

            print('---> fold=',fold)
            clf.fit(X_train, y_train)
            y_predicted= clf.predict(X_test)
            accuracy,sensitivity, specificity, precision, recall, f1, AUC=Get_model_performnace(y_test,y_predicted)
            score_clf.append(list([accuracy, sensitivity, specificity,precision, recall, f1, AUC]))

        Clf_score_clf = pd.DataFrame(np.asarray(score_clf).T)
        New_score=Clf_score_clf.mean(axis = 1, skipna = True)
        Clf_score_clf['Scores']=   list(['Accuracy','Sensitivity', 'Specificity','Precision', 'Recall','F1-score', 'ROC-AUC'])
        x = mat_filename.split("\\"); filename_sub=x[-1][:-13];filename_sub=filename_sub.replace(' ', '')
        score.append(list([name,L]+ list(New_score)))
#        clf_k.append(name'

Clf_score= pd.DataFrame(np.asarray(score).T)

Clf_score['scores']=   pd.DataFrame(list(['Classifier', 'L','Accuracy','Sensitivity', 'Specificity','Precision', 'Recall','F1-score', 'ROC-AUC']))

if save_excel==1:
    path='./Results'
    if not os.path.exists(path):
        # Create target Directory
        os.mkdir(path)

    path=path+'/'+project_name

    if not os.path.exists(path):
        # Create target Directory
        os.mkdir(path)

    Clf_score.to_csv(path+'/5Fold_CV_'+ project_name+ filename_sub+str(list_L)+'.csv', sep=',')

##%% ---------------------------- cell  ----------------------------

reason: import numpy,import scipy


text:
#!/usr/bin/env python3
import os
import time
import cv2
import pickle
import face_recognition
import numpy

# last n seconds to use for fps calc
FPS_WINDOW = 60

# if true, use face_recognition.face_distance to determin known faces
USE_FACE_DISTANCE = os.getenv('USE_FACE_DISTANCE') == '1' or False

if USE_FACE_DISTANCE:
    print('Using face_distance to determine known faces')

class FpsStats(object):
    def __init__(self):
        self.start()

    # can call start after init, or pause and start for more accuracy
    def start(self):
        self.started_at = time.time()
        self.total_frames = 0
        self.floating_frames_count = 0
        self.floating_started_at = time.time()
        self.last_floating_fps = 0

    def increment(self):
        self.total_frames += 1
        self.floating_frames_count += 1

        fps_time = time.time() - self.floating_started_at
        if fps_time > FPS_WINDOW:
            self.last_floating_fps = self.floating_frames_count / fps_time
            self.floating_started_at = time.time()
            self.floating_frames_count = 0
            print(f"fps: {self.last_floating_fps}")

    def stats(self):
        now = time.time()
        total_time = now - self.started_at
        return {
            "totalFramesRead": self.total_frames,
            "totalTime": total_time,
            "overallFps": self.total_frames / total_time,
            "fpsStartedAt": self.floating_started_at,
            "floatingFps": self.last_floating_fps
        }

print('initializing VideoCapture')
camera = cv2.VideoCapture(0)  # , apiPreference=cv2.CAP_V4L2)
if not camera.isOpened():
    raise RuntimeError('Could not start camera.')

camera.set(cv2.CAP_PROP_FRAME_WIDTH, 640)
camera.set(cv2.CAP_PROP_FRAME_HEIGHT, 480)

encodings_data = pickle.loads(
    open("data/encodings.pickle", "rb").read(), encoding='latin1')
print(f"loaded {len(encodings_data['encodings'])} encodings")

fps_stats = FpsStats()

try:
    while True:
        _, img = camera.read()
        if img is None:
            print(
                "The camera has not read data, please check whether the camera can be used normally.")
            print(
                "Use the command: 'raspistill -t 1000 -o image.jpg' to check whether the camera can be used correctly.")
            break

        fps_stats.increment()

        new_faces = face_recognition.face_locations(img)
        num_new_faces = len(new_faces)
        print(f"found {num_new_faces} faces")

        names = []
        if num_new_faces > 0:
            encodings = face_recognition.face_encodings(img, new_faces)
            for encoding in encodings:
                if USE_FACE_DISTANCE:
                    face_distances = face_recognition.face_distance(
                        encodings_data["encodings"], encoding)
                    best_match_index = numpy.argmin(face_distances)
                    if face_distances[best_match_index] < 0.65:
                        names.append(encodings_data["names"][best_match_index])
                else:
                    matches = face_recognition.compare_faces(
                        encodings_data["encodings"], encoding)
                    # check to see if we have found a match
                    if True in matches:
                        matched_indexes = [
                            i for (i, b) in enumerate(matches) if b]
                        counts = {}

                        for i in matched_indexes:
                            name = encodings_data["names"][i]
                            counts[name] = counts.get(name, 0) + 1

                        # determine the recognized face with the largest number of votes
                        names.append(max(counts, key=counts.get))

        print(f"recognized {len(names)} faces")
except KeyboardInterrupt:
    pass

print('')
print(fps_stats.stats())

ext: py

max_stars_repo_name: littlebee/shelly-bot

max_issues_repo_name: littlebee/shelly-bot

max_forks_repo_name: littlebee/shelly-bot

reason: import numpy


text:
#!/usr/bin/env python
# -*- encoding: utf-8 -*-
from __future__ import print_function

import sys
sys.path.insert(1,"../../../")  # allow us to run this standalone

import numpy as np

import h2o
from h2o.estimators.gbm import H2OGradientBoostingEstimator
from h2o.estimators.glm import H2OGeneralizedLinearEstimator
from h2o.grid import H2OGridSearch
from h2o.estimators.stackedensemble import H2OStackedEnsembleEstimator
from tests import pyunit_utils as pu

seed = 1

def test_models_not_predicting_some_classes_dont_corrupt_resulting_SE_model():

    def unique(fr):
        return np.unique(fr.as_data_frame().values).tolist()

    def scores_and_preds(models, test):
        retval = lambda d: d
        if not isinstance(models, list):
            models = [models]
            retval = lambda d: next(iter(d.values()))
        training_scores = {m.key: m.mean_per_class_error() for m in models}
        cv_scores = {m.key: m.mean_per_class_error(xval=True) for m in models}
        test_scores = {m.key: m.model_performance(test).mean_per_class_error() for m in models}
        test_predictions = {m.key: m.predict(test) for m in models}
        test_pclasses = {m.key: unique(test_predictions[m.key]['predict']) for m in models}
        return pu.ns(
            training_scores=retval(training_scores),
            cv_scores=retval(cv_scores),
            test_scores=retval(test_scores),
            test_pclasses=retval(test_pclasses),
        )

    def setup_data():
        # MNIST is multinomial classification problem
        train_full = h2o.import_file(pu.locate("bigdata/laptop/mnist/train.csv.gz"))
        # test_full = h2o.import_file(pu.locate("bigdata/laptop/mnist/test.csv.gz"))
        # train = train_full
        # test = test_full
        train, test, _ = train_full.split_frame(ratios=[.05, .1], seed=seed)
        x = train.columns[:-1]
        y = -1
        for fr in [train]:
            fr[y] = fr[y].asfactor()
        domain = unique(train[y])
        print(domain)
        return pu.ns(x=x, y=y, train=train, test=test, domain=domain)

    def train_base_models(data):
        grid = H2OGridSearch(H2OGradientBoostingEstimator,
                             search_criteria=dict(
                                 strategy='RandomDiscrete',
                                 max_models=5,
                                 seed=seed,
                             ),
                             hyper_params=dict(
                                 learn_rate=[0.5, 0.8, 1.0],
                                 max_depth=[2, 3, 4, 5],
                                 ntrees=[5, 10, 15],
                             ),
                             )
        grid.train(data.x, data.y, data.train,
                   nfolds=5,
                   fold_assignment='Modulo',
                   keep_cross_validation_predictions=True)
        return grid.models

    def train_bad_model(data):
        glm = H2OGeneralizedLinearEstimator(family='multinomial',
                                            missing_values_handling='MeanImputation',
                                            alpha=[0.0, 0.2, 0.4, 0.6, 0.8, 1.0],
                                            lambda_search=True,
                                            nfolds=5,
                                            fold_assignment='Modulo',
                                            keep_cross_validation_predictions=True,
                                            seed=seed)
        glm.train(data.x, data.y, data.train, max_runtime_secs=2)
        return glm

    def check_stackedensemble_with_AUTO_metalearner(data, models):
        se = H2OStackedEnsembleEstimator(base_models=models,
                                         metalearner_nfolds=5,
                                         seed=seed)
        se.train(data.x, data.y, data.train)
        results = scores_and_preds(se, data.test)
        print(results)
        assert data.domain == results.test_pclasses, "expected predicted classes {} but got {}".format(data.domain, results.test_pclasses)

    def check_stackedensemble_with_DRF_metalearner(data, models):
        se = H2OStackedEnsembleEstimator(base_models=models,
                                         metalearner_algorithm='DRF',
                                         metalearner_nfolds=5,
                                         seed=seed)
        se.train(data.x, data.y, data.train)
        results = scores_and_preds(se, data.test)
        print(results)
        assert data.domain == results.test_pclasses, "expected predicted classes {} but got {}".format(data.domain, results.test_pclasses)

    def check_stackedensemble_with_GLM_metalearner(data, models):
        se = H2OStackedEnsembleEstimator(base_models=models,
                                         metalearner_algorithm='GLM',
                                         metalearner_nfolds=5,
                                         seed=seed)
        se.train(data.x, data.y, data.train)
        results = scores_and_preds(se, data.test)
        print(results)
        assert data.domain != results.test_pclasses, "expected predictions not to include all target domain"
        assert len(results.test_pclasses) == 1

    def check_stackedensemble_with_GLM_metalearner_with_standardization_disabled(data, models):
        se = H2OStackedEnsembleEstimator(base_models=models,
                                         metalearner_algorithm='GLM',
                                         metalearner_nfolds=5,
                                         metalearner_params=dict(standardize=False),
                                         seed=seed)
        se.train(data.x, data.y, data.train)
        results = scores_and_preds(se, data.test)
        print(results)
        assert data.domain == results.test_pclasses, "expected predicted classes {} but got {}".format(data.domain, results.test_pclasses)

    data = setup_data()
    base_models = train_base_models(data)
    bad_model = train_bad_model(data)
    # print(scores_and_preds(bad_model, data.test))
    all_models = base_models + [bad_model]
    check_stackedensemble_with_AUTO_metalearner(data, all_models)
    check_stackedensemble_with_DRF_metalearner(data, all_models)
    check_stackedensemble_with_GLM_metalearner(data, all_models)
    check_stackedensemble_with_GLM_metalearner_with_standardization_disabled(data, all_models)

pu.run_tests([
    test_models_not_predicting_some_classes_dont_corrupt_resulting_SE_model
])


text:
"""
Tests for functions in class SolveDiffusion2D
"""

from diffusion2d import SolveDiffusion2D
import numpy as np
import pytest
import unittest

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

class TestDiffusion2D(unittest.TestCase):

    def setUp(self):
        self.solver = SolveDiffusion2D

    def test_initialize_domain(self):
        """
        Check function SolveDiffusion2D.initialize_domain
        """

        self.solver.initialize_domain(self.solver, w=100.0, h=200.0, dx=0.4, dy=0.2)

        self.assertEqual(self.solver.nx, 250)
        self.assertEqual(self.solver.ny, 1000)

    def test_initialize_physical_parameters(self):
        """
        Checks function SolveDiffusion2D.initialize_domain
        """
        self.solver.dx = 0.4
        self.solver.dy = 0.2
        self.solver.initialize_physical_parameters(self.solver, d=20.0, T_cold=500.0, T_hot=550.0)

        self.assertAlmostEqual(0.0512, self.solver.dt, 5)

    def test_set_initial_condition(self):
        """
        Checks function SolveDiffusion2D.get_initial_function
        """
        self.solver.T_cold = 500.0
        self.solver.T_hot = 550.0
        self.solver.nx = 250
        self.solver.ny = 1000
        self.solver.dx = 0.4
        self.solver.dy = 0.2

        # I feel like this is bad practice because I should calculate these
        # values manualy, but my Input values are to high, so I do it automated ;-)
        manuel_u = 500.0 * np.ones((250, 1000))
        r, cx, cy = 2, 5, 5
        r2 = r ** 2
        for i in range(250):
            for j in range(1000):
                p2 = (i * 0.4 - cx) ** 2 + (j * 0.2 - cy) ** 2
                if p2 < r2:
                    manuel_u[i, j] = 550.0

        automated_u = self.solver.set_initial_condition(self.solver)

        for i in range(self.solver.nx):
            for j in range(self.solver.ny):
                self.assertEqual(manuel_u[i,j], automated_u[i,j])


text:
import numpy as np
import pickle
from modules.pls_classifier import PLSClassifier

def labelFormat(n):
    if n == 0:
        return -1
    else:
        return n

labels = ["001", "002", "003", "004", "005", "006",
"007", "008", "009", "010"]

arq = open('LBP_files/learn.txt', 'r')

matrix_x = []
matrix_y = []

for linha in arq:
    y, x = linha.split(" | ")
    x = np.array(np.mat(x.strip("[]").strip())).ravel()
    print(y, len(x))
    matrix_x.append(x)
    matrix_y.append(y)
arq.close()

matrix_y = np.array(matrix_y)

models = []
for label in labels:
    classifier = PLSClassifier()
    y = np.char.count(matrix_y, label)
    y = list(map(labelFormat, y))
    model = classifier.fit(np.array(matrix_x), np.array(y))
    models.append(model)

path = "models_pls/lbp/"
name_arq = input("Template Name: ")
pickle.dump(models, open(path+name_arq+".sav", 'wb'))

max_stars_repo_name: AlexandreNS/butterfly-effect

max_issues_repo_name: AlexandreNS/butterfly-effect

max_forks_repo_name: AlexandreNS/butterfly-effect


text:
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from bentdna.shapefourier import ShapeAgent

class RegressAgent:

    n_begin = 0
    n_end = 9
    nlist = range(1, 10)

    def __init__(self, host, workfolder):
        self.host = host
        self.workfolder = workfolder

        self.df_an = self.__read_df_an()

        self.x_all = self.nlist_to_x(self.nlist)
        self.var_list = self.get_var_list() # y_all

        self.regr = LinearRegression()
        self.a = None

    def draw_observe_plot(self):
        fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(16,4))
        ax = axes[0]
        ax.scatter(self.nlist, self.var_list)
        ax.set_xlabel('n', fontsize=12)
        ax.set_ylabel('Var($a_n$)', fontsize=12)
        ax.set_title(self.host, fontsize=14)

        ax = axes[1]
        ax.scatter(self.x_all, self.var_list)
        ax.set_xlabel('$(1/n)^2$', fontsize=12)
        ax.set_ylabel('Var($a_n$)', fontsize=12)
        return fig, axes

    def draw_regr_var_an_nsquare_inverse(self, select_nlist):
        x = self.nlist_to_x(select_nlist)
        y = [self.get_var_an_by_n(n) for n in select_nlist]
        X_all = self.reshape_array(self.x_all)
        y_pred = self.regr.predict(X_all)

        fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8,4))
        ax.scatter(self.x_all, self.var_list,  color='black')
        ax.plot(self.x_all, y_pred, color='blue', linewidth=3)
        ax.plot(x, y, 'x', color='red')
        ax.set_xlabel('$(1/n)^2$', fontsize=12)
        ax.set_ylabel('Var($a_n$)', fontsize=12)
        return fig, ax

    def draw_regr_var_an_n(self, ax, select_nlist):
        n_continuous, y_continuous = self.get_continuous_line()
        y = [self.get_var_an_by_n(n) for n in select_nlist]

        ax.scatter(self.nlist, self.var_list)
        ax.plot(n_continuous, y_continuous, '--', color='blue', alpha=0.6)
        ax.plot(select_nlist, y, 'x', color='red')
        ax.set_xlabel('n', fontsize=12)
        ax.set_ylabel('Var($a_n$)', fontsize=12)
        ax.set_title(self.host, fontsize=14)
        return ax

    def get_continuous_line(self):
        n = np.linspace(1, 9, num=50)
        x = self.nlist_to_x(n)
        X = self.reshape_array(x)
        y = self.regr.predict(X)
        return n, y

    def do_regression(self, select_nlist):
        X, Y = self.__get_X_Y(select_nlist)
        self.regr.fit(X, Y)

        X_all = self.reshape_array(self.x_all)
        Y_all = self.reshape_array(self.var_list)
        y_pred = self.regr.predict(X_all)

        self.a = self.regr.coef_[0][0]
        mse = mean_squared_error(Y_all, y_pred)
        r2 = r2_score(Y_all, y_pred)
        print(f'Coefficients: {self.a:.3f}')
        print(f'Mean squared error: {mse:.3f}')
        print(f'Coefficient of determination: {r2:.2f}')

    def get_var_an_by_n(self, n):
        key = str(n)
        return self.df_an[key].var()

    def get_var_list(self):
        return [self.get_var_an_by_n(n) for n in self.nlist]

    def __read_df_an(self):
        s_agent = ShapeAgent(self.workfolder, self.host)
        return s_agent.read_df_an(self.n_begin, self.n_end)

    def __get_X_Y(self, select_nlist):
        x = self.nlist_to_x(select_nlist)
        y = [self.get_var_an_by_n(n) for n in select_nlist]
        X = self.reshape_array(x)
        Y = self.reshape_array(y)
        return X, Y

    def nlist_to_x(self, nlist):
        x = np.array([1/n for n in nlist])
        return np.square(x)

    def reshape_array(self, x):
        return np.matrix(x).reshape(len(x), 1)


text:
from lsd import LsdExtractor
from scipy.ndimage import gaussian_filter, convolve, maximum_filter
import h5py
import logging
import mahotas
import numpy as np
import time

logging.basicConfig(level=logging.INFO)
logging.getLogger('lsd.local_shape_descriptor').setLevel(logging.DEBUG)

k_z = np.zeros((3, 1, 1), dtype=np.float32)
k_z[0, 0, 0] = 1
k_z[2, 0, 0] = -1
k_y = np.zeros((1, 3, 1), dtype=np.float32)
k_y[0, 0, 0] = 1
k_y[0, 2, 0] = -1
k_x = np.zeros((1, 1, 3), dtype=np.float32)
k_x[0, 0, 0] = 1
k_x[0, 0, 2] = -1

k_zz = np.zeros((5, 1, 1), dtype=np.float32)
k_zz[0, 0, 0] = -1
k_zz[2, 0, 0] = 2
k_zz[4, 0, 0] = -1
k_yy = np.zeros((1, 5, 1), dtype=np.float32)
k_yy[0, 0, 0] = -1
k_yy[0, 2, 0] = 2
k_yy[0, 4, 0] = -1
k_xx = np.zeros((1, 1, 5), dtype=np.float32)
k_xx[0, 0, 0] = -1
k_xx[0, 0, 2] = 2
k_xx[0, 0, 4] = -1

k_zy = np.zeros((3, 3, 1), dtype=np.float32)
k_zy[0, 0, 0] = 1
k_zy[0, 2, 0] = -1
k_zy[2, 0, 0] = -1
k_zy[2, 2, 0] = 1
k_zx = np.zeros((3, 1, 3), dtype=np.float32)
k_zx[0, 0, 0] = 1
k_zx[0, 0, 2] = -1
k_zx[2, 0, 0] = -1
k_zx[2, 0, 2] = 1
k_yx = np.zeros((1, 3, 3), dtype=np.float32)
k_yx[0, 0, 0] = 1
k_yx[0, 0, 2] = -1
k_yx[0, 2, 0] = -1
k_yx[0, 2, 2] = 1

def create_random_segmentation(size, seed):

    np.random.seed(seed)
    peaks = np.random.random(size).astype(np.float32)
    peaks = gaussian_filter(peaks, sigma=5.0)
    max_filtered = maximum_filter(peaks, 10)
    maxima = max_filtered==peaks
    seeds, n = mahotas.label(maxima)
    print("Creating segmentation with %d segments"%n)
    return mahotas.cwatershed(1.0 - peaks, seeds).astype(np.uint64)

def get_descriptors(mask, sigma, fast=True):

    count = gaussian_filter(mask, sigma, mode='constant', output=np.float32)

    if fast:
        d_z = convolve(count, k_z, mode='constant')
        d_y = convolve(count, k_y, mode='constant')
        d_x = convolve(count, k_x, mode='constant')
        d_zz = convolve(count, k_zz, mode='constant')
        d_yy = convolve(count, k_yy, mode='constant')
        d_xx = convolve(count, k_xx, mode='constant')
        d_zy = convolve(count, k_zy, mode='constant')
        d_zx = convolve(count, k_zx, mode='constant')
        d_yx = convolve(count, k_yx, mode='constant')
    else:
        d_z = gaussian_filter(mask, sigma, order=(1, 0, 0), mode='constant', output=np.float32)
        d_y = gaussian_filter(mask, sigma, order=(0, 1, 0), mode='constant', output=np.float32)
        d_x = gaussian_filter(mask, sigma, order=(0, 0, 1), mode='constant', output=np.float32)
        d_zz = gaussian_filter(mask, sigma, order=(2, 0, 0), mode='constant', output=np.float32)
        d_yy = gaussian_filter(mask, sigma, order=(0, 2, 0), mode='constant', output=np.float32)
        d_xx = gaussian_filter(mask, sigma, order=(0, 0, 2), mode='constant', output=np.float32)
        d_zy = gaussian_filter(mask, sigma, order=(1, 1, 0), mode='constant', output=np.float32)
        d_zx = gaussian_filter(mask, sigma, order=(1, 0, 1), mode='constant', output=np.float32)
        d_yx = gaussian_filter(mask, sigma, order=(0, 1, 1), mode='constant', output=np.float32)

    d = np.stack([
        d_z, d_y, d_x,
        d_zz, d_yy, d_xx,
        d_zy, d_zx, d_yx])

    # normalize, move to [0, 1]
    count[count==0] = 1
    d = d/count + 0.5

    lsds = np.concatenate([d, count[None,:]])

    lsds[:,mask==0] = 0

    return lsds

if __name__ == "__main__":

    extractor = LsdExtractor(sigma=(5.0, 5.0, 5.0))

    segmentation = create_random_segmentation((100, 100, 100), seed=42)
    ids = np.unique(segmentation)

    start = time.time()
    lsds = extractor.get_descriptors(segmentation)
    print("Computed original LSDs in %fs"%(time.time() - start))

    start = time.time()
    lsds_compare_fast = np.zeros_like(lsds)
    for i in ids:
        lsds_compare_fast += get_descriptors(segmentation==i, (5.0, 5.0, 5.0))
    print("Computed alternative LSDs in %fs"%(time.time() - start))

    start = time.time()
    lsds_compare = np.zeros_like(lsds)
    for i in ids:
        lsds_compare += get_descriptors(segmentation==i, (5.0, 5.0, 5.0), fast=False)
    print("Computed alternative LSDs in %fs"%(time.time() - start))

    with h5py.File('test_shape_descriptor.hdf', 'w') as f:
        f['volumes/segmentation'] = segmentation
        f['volumes/lsds'] = lsds
        f['volumes/lsds_compare_fast'] = lsds_compare_fast
        f['volumes/lsds_compare'] = lsds_compare

max_stars_repo_name: pattonw/mouselight

max_issues_repo_name: pattonw/mouselight

max_forks_repo_name: pattonw/mouselight

reason: import numpy,from scipy


text:
#! /usr/bin/env python

"""
Calculating the fraction of upgoing events
CAUTION: Assuming chan0 is the uppermost and
         chan1 is below chan0
"""

import numpy as n
import pylab as p
import sys

f = open(sys.argv[1])

directions = []

for line in f.readlines():

    try:
        line = line.split()
        directions.append(float(line[3][2:-1]) - float(line[1][2:-1]))
    except:
        pass

up = 0
down = 0

for item in directions:
    if item > 0:
        down += 1
    else:
        up += 1

print "Upgoing events:", up
print "Downgoing events:", down
print "Fraction of upgoing events:", float(up)/(up+down)

max_stars_repo_name: LambdaDigamma/muonic

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: LambdaDigamma/muonic

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: LambdaDigamma/muonic

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
################################################################################
# MIT License
#
# Copyright (c) 2020 Phillip Lippe
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to conditions.
#
# Author: Deep Learning Course | Fall 2020
# Date Created: 2020-11-22
################################################################################

import torch
import torch.nn as nn
import numpy as np

class MLPEncoder(nn.Module):

    def __init__(self, input_dim=784, hidden_dims=[512], z_dim=20):
        """
        Encoder with an MLP network and ReLU activations (except the output layer).

        Inputs:
            input_dim - Number of input neurons/pixels. For MNIST, 28*28=784
            hidden_dims - List of dimensionalities of the hidden layers in the network.
                          The NN should have the same number of hidden layers as the length of the list.
            z_dim - Dimensionality of latent vector.
        """
        super().__init__()

        # For an intial architecture, you can use a sequence of linear layers and ReLU activations.
        # Feel free to experiment with the architecture yourself, but the one specified here is
        # sufficient for the assignment.
        self.z_dim = z_dim

        layers = []
        layer_sizes = [input_dim] + hidden_dims
        for idx in range(1, len(layer_sizes)):
            layers += [nn.Linear(layer_sizes[idx-1], layer_sizes[idx]),
                       nn.ELU()]
        layers += [nn.Linear(layer_sizes[-1], z_dim * 2)]
        self.encoder = nn.Sequential(*layers)

    def forward(self, x):
        """
        Inputs:
            x - Input batch with images of shape [B,C,H,W] and range 0 to 1.
        Outputs:
            mean - Tensor of shape [B,z_dim] representing the predicted mean of the latent distributions.
            log_std - Tensor of shape [B,z_dim] representing the predicted log standard deviation
                      of the latent distributions.
        """

        # Remark: Make sure to understand why we are predicting the log_std and not std

        # turn x intor [batch_size, image_vector(784)]
        x_reshape = torch.reshape(x, shape=(x.shape[0], 784))

        encoded_x = self.encoder(x_reshape)

        mean = encoded_x[:, 0:self.z_dim]
        log_std = encoded_x[:, self.z_dim:]

        return mean, log_std

class MLPDecoder(nn.Module):

    def __init__(self, z_dim=20, hidden_dims=[512], output_shape=[1, 28, 28]):
        """
        Decoder with an MLP network.
        Inputs:
            z_dim - Dimensionality of latent vector (input to the network).
            hidden_dims - List of dimensionalities of the hidden layers in the network.
                          The NN should have the same number of hidden layers as the length of the list.
            output_shape - Shape of output image. The number of output neurons of the NN must be
                           the product of the shape elements.
        """
        super().__init__()
        self.output_shape = output_shape

        output_shape = np.prod(output_shape)

        # For an intial architecture, you can use a sequence of linear layers and ReLU activations.
        # Feel free to experiment with the architecture yourself, but the one specified here is
        # sufficient for the assignment.
        layers = []
        layer_sizes = [z_dim] + hidden_dims
        for idx in range(1, len(layer_sizes)):
            layers += [nn.Linear(layer_sizes[idx-1], layer_sizes[idx]),
                       nn.ELU()]
        layers += [nn.Linear(layer_sizes[-1], output_shape)]
        self.decoder = nn.Sequential(*layers)

    def forward(self, z):
        """
        Inputs:
            z - Latent vector of shape [B,z_dim]
        Outputs:
            x - Prediction of the reconstructed image based on z.
                This should be a logit output *without* a sigmoid applied on it.
                Shape: [B,output_shape[0],output_shape[1],output_shape[2]]
        """
        # get batch_size for reshape
        batch_size = z.shape[0]

        # final output shape
        output_shape = tuple([batch_size] + self.output_shape)

        # decode z
        decoded_x = self.decoder(z)

        # reshape
        x = torch.reshape(decoded_x, shape=output_shape)

        return x

    @property
    def device(self):
        """
        Property function to get the device on which the decoder is.
        Might be helpful in other functions.
        """
        return next(self.parameters()).device


text:
import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy

from ..measure import ConditionedLognormalSampler

class ScalarImage:
    """
    Class containing a scalar image.
    """

    def __init__(self, height=1000, width=1000):
        """ Instantiate scalar image with shape (<height>, <width>). """
        self.height = height
        self.width = width
        self.initialize()

    @property
    def shape(self):
        """ Image shape. """
        return self.im.shape[-2:]

    @property
    def pixels(self):
        """ Returns image pixels. """
        return self.im.ravel()

    @property
    def max(self):
        """ Maximum pixel intensity. """
        return self.im.max()

    @property
    def im_normalized(self):
        """ Image normalized by the maximum value. """
        return self.im/self.max

    def percentile(self, q):
        """ 98th percentile of pixel intensities. """
        return np.percentile(self.im.ravel(), q=q)

    def initialize(self):
        """ Initialize blank image. """
        self.im = np.zeros((self.height, self.width), dtype=np.float64)

    def fill(self, mu=0.1, sigma=0.1):
        """
        Fill image background with values sampled from a lognormal distribution.

        Args:

            mu (float) - mean of underlying normal distribution

            sigma (float) - std dev of underlying normal distribution

        """
        pixels = np.exp(np.random.normal(np.log(mu), sigma, size=self.shape))
        self.im[:, :] = pixels

    @staticmethod
    def _render(im, vmin=0, vmax=None, cmap=plt.cm.Greys, size=5, ax=None):
        """
        Render image.

        Args:

            im (np.ndarray[float]) - image

            vmin, vmax (int) - colormap bounds

            cmap (matplotlib.ColorMap or str) - if value is 'r', 'g', or 'b', use RGB colorscheme

            size (int) - image panel size, in inches

            ax (matplotlib.axes.AxesSubplot) - if None, create figure

        """

        if ax is None:
            fig, ax = plt.subplots(figsize=(size, size))

        if vmax is None:
            vmax = im.max()

        # render image
        if type(cmap) == str:
            assert cmap in 'rgb', 'Color not recognized.'
            im_rgb = np.zeros(im.shape+(3,), dtype=np.float64)
            im_rgb[:,:,'rgb'.index(cmap)] = (im-vmin)/(vmax-vmin)
            im_rgb[im_rgb>1.] = 1.
            ax.imshow(im_rgb)
        else:
            ax.imshow(im, vmin=vmin, vmax=vmax, cmap=cmap)

        # invert axis and remove ticks
        ax.invert_yaxis()
        ax.axis('off')

    def render(self, **kwargs):
        """ Render image. """
        self._render(self.im.T, **kwargs)

    def render_blank(self, **kwargs):
        """ Render image. """
        self._render(np.zeros(self.shape, dtype=int), **kwargs)

    def center_xycoords(self, xy, shrinkage=0.9):
        """ Project zero-centered coordinates to center of image. """
        center_x, center_y = self.shape[0]/2,  self.shape[1]/2
        centered_xy = deepcopy(xy)
        centered_xy[:, 0] = ((xy[:, 0]*center_x*shrinkage) + center_x)
        centered_xy[:, 1] = ((xy[:, 1]*center_y*shrinkage) + center_y)
        return centered_xy.astype(int)

class DependentScalarImage(ScalarImage):
    """
    Class defines a scalar image whose pixel intensities are sampled with some dependence upon another scalar image.
    """

    def __init__(self, pixels, mean, sigma):
        """ Instantiate a dependent scalar image. """
        super().__init__(*pixels.shape)
        x = np.log(pixels.ravel())
        self.sampler = ConditionedLognormalSampler(x, np.log(mean), sigma)

    def fill(self, rho=0.0):
        """ Generate randomly sampled pixel values. """
        pixels = self.sampler.sample(rho=rho)
        self.im[:, :] = pixels.reshape(self.shape)

max_stars_repo_name: sebastianbernasek/growth

max_issues_repo_name: sbernasek/growth

max_forks_repo_name: sbernasek/growth


text:
import sys, math
import numpy as np

class Point3D:
    def __init__(self, x = 0, y = 0, z = 0):
        self.x, self.y, self.z = float(x), float(y), float(z)

    def view_calc(self, camera_position):
        #CAMERA CO-ORDINATES

        eye = np.array(camera_position)
        target = np.array([0, 0, 0])
        up = np.array([0, 1, 0])
        if eye[0] == 0 and eye[2] == 0:
            up = np.array([0, 0, -1])

        vz = eye - target
        vz = vz / np.linalg.norm(vz)
        vx = np.cross(up, vz)
        vx = vx / np.linalg.norm(vx)
        vy = np.cross(vz, vx)

        view = np.matrix([[vx[0], vx[1], vx[2], 0], \
                          [vy[0], vy[1], vy[2], 0], \
                          [vz[0], vz[1], vz[2], 0], \
                          [-np.dot(vx,eye), -np.dot(vy,eye), -np.dot(vz,eye), 1]])

        return view

    def perspective2(self, win_width, win_height):

        near = 0.001
        far = 100.0
        angle_of_view = 45.0
        aspect_ratio = win_width / win_height

        size = near * math.tan(math.radians(angle_of_view) / 2.0)
        left = -size
        right = size
        bottom = -size / aspect_ratio
        top = size / aspect_ratio

        proj = np.matrix([[2*near/(right - left), 0.0, (right+left)/(right-left), 0.0], \
                          [0.0, 2*near/(top - bottom), (top+bottom)/(top-bottom), 0.0], \
                          [0.0, 0.0, -1.0*(far+near)/(far-near), (-2.0*far*near)/(far-near)], \
                          [0.0, 0.0, -1.0, 0.0]])
        return proj

    def perspective(self, x, y, z, win_width, win_height):

        fov = 512
        viewer_distance = 512
        factor = fov / (viewer_distance + z)
        x = -x * factor + win_width / 2
        y = -y * factor + win_height / 2
        return [x, y, 1]

    def project(self, win_width, win_height, camera_position):
        """ Transforms this 3D point to 2D using a perspective projection. """

        # START CALCULATE SCREEN COORDINATES
        point = np.matrix([[self.x], \
                           [self.y], \
                           [self.z], \
                           [1.]])

        view = self.view_calc(camera_position)

        # Scale Part of Model Matrix
        scale = 1
        scalemat = np.matrix([[scale, 0, 0, 0], \
                             [0, scale, 0, 0], \
                             [0, 0, scale, 0], \
                             [0, 0, 0, 1]])

        self.modelView = view * scalemat
        self.eyeView = self.modelView * point
        final = self.perspective(self.eyeView[0,0], self.eyeView[1,0], self.eyeView[2,0], win_width, win_height)
        scrn_x = final[0]
        scrn_y = final[1]

        x, y = scrn_x, scrn_y
        return Point3D(x, y, 1)


text:
import numpy as np

import pandas as pd

from .utils import raises
from .. import ValidationError
from ..validators import broadcastable, check_type, dict_validator, has, \
    nonsingular, true

class A(object):
    pass

def test_check_type():
    func = check_type(A)
    a = A()
    func(a)
    assert raises(ValidationError, func, 1)
    func = check_type((dict, list))
    assert raises(ValidationError, func, (1, 2))
    func([1, 2])
    func({'a': 1})

def test_nonsingular():
    assert raises(ValidationError, nonsingular, [[0, 0], [0, 0]])
    nonsingular([[0, 1], [1, 0]])

def test_broadcastable():
    data = {'a': np.random.random((2, 1)),
            'b': np.random.random((3, 1)),
            'c': 3}
    func = broadcastable('a', 'b')
    assert raises(ValidationError, func, data)
    func = broadcastable('a', 'c')
    func(data)

def test_has():
    data = pd.DataFrame({'a': [1, 2, 3, 4, 5],
                         'b': [1, 2, 3, 4, 5],
                         'c': [1, 2, 3, 4, 5]})
    func = has('a', 'b')
    func(data)
    func = has('a', 'dd')
    assert raises(ValidationError, func, data)

def test_dict_validator():
    good_dict = {'a': 1, 'b': 2, 'c': [1, 2, 3, 4]}
    bad_dict = {'a': 1, 'c': [1, 2, 3, 4]}
    validators = dict(b=lambda x: true(x is not None, "must be present"))
    func = dict_validator(validators)
    func(good_dict)
    assert raises(ValidationError, func, bad_dict)

max_stars_repo_name: hhuuggoo/thedoctor

max_issues_repo_name: hhuuggoo/thedoctor

max_forks_repo_name: hhuuggoo/thedoctor


text:
from PyQt5 import QtWidgets
from PyQt5.QtGui import QPainter
from PyQt5.QtWidgets import QStyle, QStyleOptionSlider
from PyQt5.QtCore import QRect, QPoint, Qt
from sklearn.datasets import make_blobs, make_moons, make_circles
import time
from matplotlib.pyplot import Polygon
from scipy.spatial import ConvexHull
from matplotlib import colors
import numpy as np

def choose_dataset(chosen_dataset, n_points):
    X = None

    if chosen_dataset == "blobs":
        X = make_blobs(
            n_samples=n_points,
            centers=4,
            n_features=2,
            cluster_std=1.5,
            random_state=42,
        )[0]
    elif chosen_dataset == "moons":
        X = make_moons(n_samples=n_points, noise=0.05, random_state=42)[0]
    elif chosen_dataset == "scatter":
        X = make_blobs(
            n_samples=n_points, cluster_std=[2.5, 2.5, 2.5], random_state=42
        )[0]
    elif chosen_dataset == "circle":
        X = make_circles(n_samples=n_points, noise=0, random_state=42)[0]

    return X

def pause_execution(seconds):
    time.sleep(seconds)

def encircle(x, y, ax, **kw):
    p = np.c_[x, y]
    hull = ConvexHull(p)
    poly = Polygon(p[hull.vertices, :], **kw)
    ax.add_patch(poly)

def convert_colors(dict_colors, alpha=0.5):
    new_dict_colors = {}

    for i, col in enumerate(dict_colors.values()):
        new_dict_colors[i] = tuple(list(colors.to_rgb(col)) + [alpha])

    return new_dict_colors

class LabeledSlider(QtWidgets.QWidget):
    def __init__(
        self,
        minimum,
        maximum,
        interval=1,
        single_step=1,
        orientation=Qt.Horizontal,
        labels=None,
        parent=None,
    ):
        super(LabeledSlider, self).__init__(parent=parent)

        levels = range(minimum, maximum + 1, interval)
        if labels is not None:
            if not isinstance(labels, (tuple, list)):
                raise Exception("<labels> is a list or tuple.")
            if len(labels) != len(levels):
                raise Exception("Size of <labels> doesn't match levels.")
            self.levels = list(zip(levels, labels))
        else:
            self.levels = list(zip(levels, map(str, levels)))

        if orientation == Qt.Horizontal:
            self.layout = QtWidgets.QVBoxLayout(self)
        elif orientation == Qt.Vertical:
            self.layout = QtWidgets.QHBoxLayout(self)
        else:
            raise Exception("<orientation> wrong.")

        # gives some space to print labels
        self.left_margin = 0
        self.top_margin = 0
        self.right_margin = 0
        self.bottom_margin = 0

        self.layout.setContentsMargins(
            self.left_margin, self.top_margin, self.right_margin, self.bottom_margin
        )

        self.sl = QtWidgets.QSlider(orientation, self)
        self.sl.setMinimum(minimum)
        self.sl.setMaximum(maximum)
        self.sl.setValue(minimum)
        self.sl.setSingleStep(single_step)
        if orientation == Qt.Horizontal:
            self.sl.setTickPosition(QtWidgets.QSlider.TicksBelow)
            self.sl.setMinimumWidth(10)  # just to make it easier to read
        else:
            self.sl.setTickPosition(QtWidgets.QSlider.TicksLeft)
            self.sl.setMinimumHeight(300)  # just to make it easier to read
        self.sl.setTickInterval(interval)

        self.layout.addWidget(self.sl)

    def paintEvent(self, e):

        super(LabeledSlider, self).paintEvent(e)

        style = self.sl.style()
        painter = QPainter(self)
        st_slider = QStyleOptionSlider()
        st_slider.initFrom(self.sl)
        st_slider.orientation = self.sl.orientation()

        length = style.pixelMetric(QStyle.PM_SliderLength, st_slider, self.sl)
        available = style.pixelMetric(
            QStyle.PM_SliderSpaceAvailable, st_slider, self.sl
        )

        for v, v_str in self.levels:

            # get the size of the label
            rect = painter.drawText(QRect(), Qt.TextDontPrint, v_str)

            if self.sl.orientation() == Qt.Horizontal:
                # I assume the offset is half the length of slider, therefore
                # + length//2
                x_loc = (
                    QStyle.sliderPositionFromValue(
                        self.sl.minimum(), self.sl.maximum(), v, available
                    )
                    + length // 2
                )

                # left bound of the text = center - half of text width + L_margin
                left = x_loc - rect.width() // 2 + self.left_margin
                bottom = self.rect().bottom()

                # enlarge margins if clipping
                # if v == self.sl.minimum():
                #     if left <= 0:
                #         self.left_margin = rect.width() // 2 - x_loc
                #     if self.bottom_margin <= rect.height():
                #         self.bottom_margin = rect.height()
                #
                #     self.layout.setContentsMargins(self.left_margin,
                #                                    self.top_margin, self.right_margin,
                #                                    self.bottom_margin)
                #
                # if v == self.sl.maximum() and rect.width() // 2 >= self.right_margin:
                #     self.right_margin = rect.width() // 2
                #     self.layout.setContentsMargins(self.left_margin,
                #                                    self.top_margin, self.right_margin,
                #                                    self.bottom_margin)

            else:
                y_loc = QStyle.sliderPositionFromValue(
                    self.sl.minimum(), self.sl.maximum(), v, available, upsideDown=True
                )

                bottom = y_loc + length // 2 + rect.height() // 2 + self.top_margin - 3
                # there is a 3 px offset that I can't attribute to any metric

                left = self.left_margin - rect.width()
                if left <= 0:
                    self.left_margin = rect.width() + 2
                    self.layout.setContentsMargins(
                        self.left_margin,
                        self.top_margin,
                        self.right_margin,
                        self.bottom_margin,
                    )

            pos = QPoint(left, bottom)
            painter.drawText(pos, v_str)

        return

max_stars_repo_name: guglielmosanchini/ClustVizGUI

max_issues_repo_name: guglielmosanchini/ClustVizGUI

max_forks_repo_name: guglielmosanchini/ClustVizGUI


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

from pySDC.helpers.stats_helper import filter_stats, sort_stats
from pySDC.implementations.collocation_classes.gauss_lobatto import CollGaussLobatto
from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order
from pySDC.playgrounds.Boris.spiraling_particle_HookClass import particles_output
from pySDC.playgrounds.Boris.spiraling_particle_ProblemClass import planewave_single

def main(dt, Tend):

    # initialize level parameters
    level_params = dict()
    level_params['restol'] = 1E-08
    level_params['dt'] = dt

    # initialize sweeper parameters
    sweeper_params = dict()
    sweeper_params['collocation_class'] = CollGaussLobatto
    sweeper_params['num_nodes'] = 3

    # initialize problem parameters for the Penning trap
    problem_params = dict()
    problem_params['delta'] = 1
    problem_params['a0'] = 0.07
    problem_params['u0'] = np.array([[0, -1, 0], [0.05, 0.01, 0], [1], [1]])

    # initialize step parameters
    step_params = dict()
    step_params['maxiter'] = 20

    # initialize controller parameters
    controller_params = dict()
    controller_params['hook_class'] = particles_output  # specialized hook class for more statistics and output
    controller_params['logger_level'] = 30
    # controller_params['log_to_file'] = True
    # controller_params['fname'] = 'step_3_B_out.txt'

    # Fill description dictionary for easy hierarchy creation
    description = dict()
    description['problem_class'] = planewave_single
    description['problem_params'] = problem_params
    description['sweeper_class'] = boris_2nd_order
    description['sweeper_params'] = sweeper_params
    description['level_params'] = level_params
    # description['space_transfer_class'] = particles_to_particles # this is only needed for more than 2 levels
    description['step_params'] = step_params

    # instantiate the controller (no controller parameters used here)
    controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)

    # set time parameters
    t0 = 0.0
    Tend = Tend

    # get initial values on finest level
    P = controller.MS[0].levels[0].prob
    uinit = P.u_init()

    # call main function to get things done...
    uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)

    return uinit, stats, problem_params['a0']

def plot_error_and_positions(uinit, stats, a0):

    extract_stats = filter_stats(stats, type='energy')
    sortedlist_stats = sort_stats(extract_stats, sortby='time')

    R0 = np.linalg.norm(uinit.pos[:])
    H0 = 1 / 2 * np.dot(uinit.vel[:].T, uinit.vel[:]) + a0 / R0

    energy_err = [abs(entry[1] - H0) / H0 for entry in sortedlist_stats]

    plt.figure()
    plt.plot(energy_err, 'bo--')

    plt.xlabel('Time')
    plt.ylabel('Error in hamiltonian')

    plt.savefig('spiraling_particle_error_ham.png',  transparent=True, bbox_inches='tight')

    extract_stats = filter_stats(stats, type='position')
    sortedlist_stats = sort_stats(extract_stats, sortby='time')

    xpositions = [item[1][0] for item in sortedlist_stats]
    ypositions = [item[1][1] for item in sortedlist_stats]

    plt.figure()
    plt.xlim([-1.5, 1.5])
    plt.ylim([-1.5, 1.5])
    plt.xlabel('x')
    plt.ylabel('y')

    plt.scatter(xpositions, ypositions)
    plt.savefig('spiraling_particle_positons.png',  transparent=True, bbox_inches='tight')

    plt.show()

if __name__ == "__main__":
    uinit, stats, a0 = main(dt=1.0, Tend=5000)
    plot_error_and_positions(uinit, stats, a0)

max_stars_repo_name: brownbaerchen/pySDC

max_issues_repo_name: brownbaerchen/pySDC

max_forks_repo_name: brownbaerchen/pySDC


text:
import numpy as np
from qutip import basis
from qutip_qip.circuit import QubitCircuit
from qutip_qip.device import LinearSpinChain, SpinChainModel, Processor, ModelProcessor
from qutip_qip.compiler import SpinChainCompiler
from qutip_qip.noise import RelaxationNoise
qc = QubitCircuit(3)
qc.add_gate("X", targets=2)
qc.add_gate("SNOT", targets=0)
qc.add_gate("SNOT", targets=1)
qc.add_gate("SNOT", targets=2)

# Oracle function f(x)
qc.add_gate(
    "CNOT", controls=0, targets=2)
qc.add_gate(
    "CNOT", controls=1, targets=2)

qc.add_gate("SNOT", targets=0)
qc.add_gate("SNOT", targets=1)

init_state = basis([2,2,2], [0,0,0])
final_state = qc.run(init_state)

processor = LinearSpinChain(
    num_qubits=3, sx=0.25)

processor.load_circuit(qc)

tlist = np.linspace(0, 20, 300)
result = processor.run_state(
    init_state, tlist=tlist)

processor.add_noise(
    RelaxationNoise(t2=30))

# Section 4.1 Model
model = SpinChainModel(
    num_qubits=3, setup="circular", g=1)
processor = Processor(model=model)

model.get_control(label="sx0")

model.get_control_labels()

# Section 4.3 Scheduler
from qutip_qip.compiler import Scheduler, Instruction

Scheduler("ASAP").schedule(qc)

inst_list = []
for gate in qc.gates:
    if gate.name in ("SNOT", "X"):
        inst_list.append(
            Instruction(gate, duration=1
            )
        )
    else:
        inst_list.append(
            Instruction(gate, duration=2
            )
        )
Scheduler("ALAP").schedule(inst_list)

max_stars_repo_name: bopardikarsoham/qutip-qip

max_issues_repo_name: bopardikarsoham/qutip-qip

max_forks_repo_name: BoxiLi/qutip-qip


text:
# Copyright 2019 NREL

# Licensed under the Apache License, Version 2.0 (the "License"); you may not use
# this file except in compliance with the License. You may obtain a copy of the
# License at http://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software distributed
# under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR
# CONDITIONS OF ANY KIND, either express or implied. See the License for the
# specific language governing permissions and limitations under the License.

import numpy as np
import sys
import datetime
from wisdem.ccblade import CCAirfoil, CCBlade
from scipy import interpolate, gradient, integrate

# Some useful constants
now = datetime.datetime.now()
pi = np.pi
rad2deg = np.rad2deg(1)
deg2rad = np.deg2rad(1)
rpm2RadSec = 2.0*(np.pi)/60.0
RadSec2rpm = 60/(2.0 * np.pi)

class Controller():
    """
    Class Controller used to calculate controller tunings parameters

    Methods:
    -------
    tune_controller

    Parameters:
    -----------
    controller_params: dict
                       Dictionary containing controller paramaters that need to be defined
    """

    def __init__(self, controller_params):
        '''
        Load controller tuning parameters from input dictionary
        '''

        print('-----------------------------------------------------------------------------')
        print('   Tuning a reference wind turbine controller using NREL\'s ROSCO toolbox    ')
        # print('      Developed by Nikhar J. Abbas for collaborative research purposes.      ')
        print('-----------------------------------------------------------------------------')

        # Controller Flags
        self.LoggingLevel = controller_params['LoggingLevel']
        self.F_LPFType = controller_params['F_LPFType']
        self.F_NotchType = controller_params['F_NotchType']
        self.IPC_ControlMode = controller_params['IPC_ControlMode']
        self.VS_ControlMode = controller_params['VS_ControlMode']
        self.PC_ControlMode = controller_params['PC_ControlMode']
        self.Y_ControlMode = controller_params['Y_ControlMode']
        self.SS_Mode = controller_params['SS_Mode']
        self.WE_Mode = controller_params['WE_Mode']
        self.PS_Mode = controller_params['PS_Mode']
        self.SD_Mode = controller_params['SD_Mode']
        self.Fl_Mode = controller_params['Fl_Mode']
        self.Flp_Mode = controller_params['Flp_Mode']

        # Necessary parameters
        self.zeta_pc = controller_params['zeta_pc']
        self.omega_pc = controller_params['omega_pc']
        self.zeta_vs = controller_params['zeta_vs']
        self.omega_vs = controller_params['omega_vs']
        if self.Flp_Mode > 0:
            self.zeta_flp = controller_params['zeta_flp']
            self.omega_flp = controller_params['omega_flp']

        # Optional parameters, default to standard if not defined
        if controller_params['min_pitch']:
            self.min_pitch = controller_params['min_pitch']
        else:
            self.min_pitch = None

        if controller_params['max_pitch']:
            self.max_pitch = controller_params['max_pitch']
        else:
            self.max_pitch = 90*deg2rad      # Default to 90 degrees max pitch

        if controller_params['vs_minspd']:
            self.vs_minspd = controller_params['vs_minspd']
        else:
            self.vs_minspd = None

        if controller_params['ss_vsgain']:
            self.ss_vsgain = controller_params['ss_vsgain']
        else:
            self.ss_vsgain = 1.      # Default to 100% setpoint shift

        if controller_params['ss_pcgain']:
            self.ss_pcgain = controller_params['ss_pcgain']
        else:
            self.ss_pcgain = 0.001      # Default to 0.1% setpoint shift

        if controller_params['ss_cornerfreq']:
            self.ss_cornerfreq = controller_params['ss_cornerfreq']
        else:
            self.ss_cornerfreq = .62831850001     # Default to 10 second time constant

        if controller_params['ps_percent']:
            self.ps_percent = controller_params['ps_percent']
        else:
            self.ps_percent = 0.75      # Default to 75% peak shaving

        # critical damping if LPFType = 2
        if controller_params['F_LPFType']:
            if controller_params['F_LPFType'] == 2:
                self.F_LPFDamping = 0.7
            else:
                self.F_LPFDamping = 0.0
        else:
            self.F_LPFDamping = 0.0

        # Shutdown filter default cornering freq at 15s time constant
        if controller_params['sd_cornerfreq']:
            self.sd_cornerfreq = controller_params['sd_cornerfreq']
        else:
            self.sd_cornerfreq = 0.41888

        if controller_params['sd_maxpit']:
            self.sd_maxpit = controller_params['sd_maxpit']
        else:
            self.sd_maxpit = None

        if controller_params['flp_maxpit']:
            self.flp_maxpit = controller_params['flp_maxpit']
        else:
            if controller_params['Flp_Mode'] > 0:
                self.flp_maxpit = 10.0 * deg2rad
            else:
                self.flp_maxpit = 0.0

    def tune_controller(self, turbine):
        """
        Given a turbine model, tune a controller based on the NREL generic controller tuning process

        Parameters:
        -----------
        turbine : class
                  Turbine class containing necessary turbine information to accurately tune the controller.
        """
        # -------------Load Parameters ------------- #
        # Re-define Turbine Parameters for shorthand
        J = turbine.J                           # Total rotor inertial (kg-m^2)
        rho = turbine.rho                       # Air density (kg/m^3)
        R = turbine.rotor_radius                    # Rotor radius (m)
        Ar = np.pi*R**2                         # Rotor area (m^2)
        Ng = turbine.Ng                         # Gearbox ratio (-)
        rated_rotor_speed = turbine.rated_rotor_speed               # Rated rotor speed (rad/s)

        # -------------Define Operation Points ------------- #
        TSR_rated = rated_rotor_speed*R/turbine.v_rated  # TSR at rated

        # separate wind speeds by operation regions
        v_below_rated = np.arange(turbine.v_min,turbine.v_rated,0.5)             # below rated
        v_above_rated = np.arange(turbine.v_rated+0.5,turbine.v_max,0.5)             # above rated
        v = np.concatenate((v_below_rated, v_above_rated))

        # separate TSRs by operations regions
        TSR_below_rated = np.ones(len(v_below_rated))*turbine.TSR_operational # below rated
        TSR_above_rated = rated_rotor_speed*R/v_above_rated                     # above rated
        TSR_op = np.concatenate((TSR_below_rated, TSR_above_rated))   # operational TSRs

        # Find expected operational Cp values
        Cp_above_rated = turbine.Cp.interp_surface(0,TSR_above_rated[0])             # Cp during rated operation (not optimal). Assumes cut-in bld pitch to be 0
        Cp_op_br = np.ones(len(v_below_rated)) * turbine.Cp.max              # below rated
        Cp_op_ar = Cp_above_rated * (TSR_above_rated/TSR_rated)**3           # above rated
        Cp_op = np.concatenate((Cp_op_br, Cp_op_ar))                # operational CPs to linearize around
        pitch_initial_rad = turbine.pitch_initial_rad
        TSR_initial = turbine.TSR_initial

        # initialize variables
        pitch_op    = np.empty(len(TSR_op))
        dCp_beta    = np.empty(len(TSR_op))
        dCp_TSR     = np.empty(len(TSR_op))
        dCt_beta    = np.empty(len(TSR_op))
        dCt_TSR     = np.empty(len(TSR_op))
        Ct_op       = np.empty(len(TSR_op))

        # ------------- Find Linearized State "Matrices" ------------- #
        for i in range(len(TSR_op)):
            # Find pitch angle as a function of expected operating CP for each TSR
            Cp_TSR = np.ndarray.flatten(turbine.Cp.interp_surface(turbine.pitch_initial_rad, TSR_op[i]))     # all Cp values for a given tsr
            Cp_op[i] = np.clip(Cp_op[i], np.min(Cp_TSR), np.max(Cp_TSR))        # saturate Cp values to be on Cp surface
            f_cp_pitch = interpolate.interp1d(Cp_TSR,pitch_initial_rad)         # interpolate function for Cp(tsr) values
            pitch_op[i] = f_cp_pitch(Cp_op[i])                                  # expected operation blade pitch values
            dCp_beta[i], dCp_TSR[i] = turbine.Cp.interp_gradient(pitch_op[i],TSR_op[i])       # gradients of Cp surface in Beta and TSR directions
            dCt_beta[i], dCt_TSR[i] = turbine.Ct.interp_gradient(pitch_op[i],TSR_op[i])       # gradients of Cp surface in Beta and TSR directions

            # Thrust
            Ct_TSR      = np.ndarray.flatten(turbine.Ct.interp_surface(turbine.pitch_initial_rad, TSR_op[i]))     # all Cp values for a given tsr
            f_ct        = interpolate.interp1d(pitch_initial_rad,Ct_TSR)
            Ct_op[i]    = f_ct(pitch_op[i])
            Ct_op[i]    = np.clip(Ct_op[i], np.min(Ct_TSR), np.max(Ct_TSR))        # saturate Ct values to be on Ct surface

        # Define minimum pitch saturation to be at Cp-maximizing pitch angle if not specifically defined
        if not self.min_pitch:
            self.min_pitch = 0.0
            self.min_pitch = max(self.min_pitch,pitch_op[0])

        # Full Cx surface gradients
        dCp_dbeta   = dCp_beta/np.diff(pitch_initial_rad)[0]
        dCp_dTSR    = dCp_TSR/np.diff(TSR_initial)[0]
        dCt_dbeta   = dCt_beta/np.diff(pitch_initial_rad)[0]
        dCt_dTSR    = dCt_TSR/np.diff(TSR_initial)[0]

        # Linearized system derivatives
        dtau_dbeta      = Ng/2*rho*Ar*R*(1/TSR_op)*dCp_dbeta*v**2
        dtau_dlambda    = Ng/2*rho*Ar*R*v**2*(1/(TSR_op**2))*(dCp_dTSR*TSR_op - Cp_op)
        dlambda_domega  = R/v/Ng
        dtau_domega     = dtau_dlambda*dlambda_domega

        dlambda_dv      = -(TSR_op/v)

        Pi_beta         = 1/2 * rho * Ar * v**2 * dCt_dbeta
        Pi_omega        = 1/2 * rho * Ar * R * v * dCt_dTSR
        Pi_wind         = 1/2 * rho * Ar * v**2 * dCt_dTSR * dlambda_dv + rho * Ar * v * Ct_op

        # Second order system coefficients
        A = dtau_domega/J             # Plant pole
        B_tau = -Ng**2/J              # Torque input
        B_beta = dtau_dbeta/J         # Blade pitch input

        # Wind Disturbance Input
        dtau_dv = (0.5 * rho * Ar * 1/rated_rotor_speed) * (dCp_dTSR*dlambda_dv*v**3 + Cp_op*3*v**2)
        B_wind = dtau_dv/J # wind speed input - currently unused

        # separate and define below and above rated parameters
        A_vs = A[0:len(v_below_rated)]          # below rated
        A_pc = A[len(v_below_rated):len(v)]     # above rated
        B_tau = B_tau * np.ones(len(v))

        # -- Find gain schedule --
        self.pc_gain_schedule = ControllerTypes()
        self.pc_gain_schedule.second_order_PI(self.zeta_pc, self.omega_pc,A_pc,B_beta[len(v_below_rated):len(v)],linearize=True,v=v_above_rated)
        self.vs_gain_schedule = ControllerTypes()
        self.vs_gain_schedule.second_order_PI(self.zeta_vs, self.omega_vs,A_vs,B_tau[0:len(v_below_rated)],linearize=False,v=v_below_rated)

        # -- Find K for Komega_g^2 --
        self.vs_rgn2K = (pi*rho*R**5.0 * turbine.Cp.max) / (2.0 * turbine.Cp.TSR_opt**3 * Ng**3)
        self.vs_refspd = min(turbine.TSR_operational * turbine.v_rated/R, turbine.rated_rotor_speed) * Ng

        # -- Define some setpoints --
        # minimum rotor speed saturation limits
        if self.vs_minspd:
            self.vs_minspd = np.maximum(self.vs_minspd, (turbine.TSR_operational * turbine.v_min / turbine.rotor_radius) * Ng)
        else:
            self.vs_minspd = (turbine.TSR_operational * turbine.v_min / turbine.rotor_radius) * Ng
        self.pc_minspd = self.vs_minspd

        # max pitch angle for shutdown
        if self.sd_maxpit:
            self.sd_maxpit = self.sd_maxpit
        else:
            self.sd_maxpit = pitch_op[-1]

        # Store some variables
        self.v              = v                                  # Wind speed (m/s)
        self.v_below_rated  = v_below_rated
        self.pitch_op       = pitch_op
        self.pitch_op_pc    = pitch_op[len(v_below_rated):len(v)]
        self.TSR_op         = TSR_op
        self.A              = A
        self.B_beta         = B_beta
        self.B_tau          = B_tau
        self.B_wind         = B_wind
        self.TSR_op         = TSR_op
        self.omega_op       = TSR_op*v/R
        self.Pi_omega       = Pi_omega
        self.Pi_beta        = Pi_beta
        self.Pi_wind        = Pi_wind

        # - Might want these to debug -
        # self.Cp_op = Cp_op

        # --- Minimum pitch saturation ---
        self.ps_min_bld_pitch = np.ones(len(self.pitch_op)) * self.min_pitch
        self.ps = ControllerBlocks()

        if self.PS_Mode == 1:  # Peak Shaving
            self.ps.peak_shaving(self, turbine)
        elif self.PS_Mode == 2: # Cp-maximizing minimum pitch saturation
            self.ps.min_pitch_saturation(self,turbine)
        elif self.PS_Mode == 3: # Peak shaving and Cp-maximizing minimum pitch saturation
            self.ps.peak_shaving(self, turbine)
            self.ps.min_pitch_saturation(self,turbine)

        # --- Floating feedback term ---
        if self.Fl_Mode == 1: # Floating feedback
            Kp_float = (dtau_dv/dtau_dbeta) * turbine.TowerHt * Ng
            self.Kp_float = Kp_float[len(v_below_rated)]
            # Turn on the notch filter if floating
            self.F_NotchType = 2

            # And check for .yaml input inconsistencies
            if turbine.twr_freq == 0.0 or turbine.ptfm_freq == 0.0:
                print('WARNING: twr_freq and ptfm_freq should be defined for floating turbine control!!')
        else:
            self.Kp_float = 0.0

        # Flap actuation
        if self.Flp_Mode >= 1:
            self.flp_angle = 0.0
            try:
                self.tune_flap_controller(turbine)
            except AttributeError:
                print('ERROR: If Flp_Mode > 0, you need to have blade information loaded in the turbine object.')
                raise
            except UnboundLocalError:
                print('ERROR: You are attempting to tune a flap controller for a blade without flaps!')
                raise
        else:
            self.flp_angle = 0.0
            self.Ki_flap = np.array([0.0])
            self.Kp_flap = np.array([0.0])

    def tune_flap_controller(self,turbine):
        '''
        Tune controller for distributed aerodynamic control

        Parameters:
        -----------
        turbine : class
                  Turbine class containing necessary turbine information to accurately tune the controller.
        '''
        # Find blade aerodynamic coefficients
        v_rel = []
        phi_vec = []
        alpha=[]
        for i, _ in enumerate(self.v):
            turbine.cc_rotor.induction_inflow=True
            # Axial and tangential inductions
            a, ap, alpha0, cl, cd = turbine.cc_rotor.distributedAeroLoads(self.v[i], self.omega_op[i], self.pitch_op[i], 0.0)
            # Relative windspeed
            v_rel.append([np.sqrt(self.v[i]**2*(1-a)**2 + self.omega_op[i]**2*turbine.span**2*(1-ap)**2)])
            # Inflow wind direction
            phi_vec.append(self.pitch_op[i] + turbine.twist*deg2rad)

        # Lift and drag coefficients
        Cl0 = np.zeros_like(turbine.af_data)
        Cd0 = np.zeros_like(turbine.af_data)
        Clp = np.zeros_like(turbine.af_data)
        Cdp = np.zeros_like(turbine.af_data)
        Clm = np.zeros_like(turbine.af_data)
        Cdm = np.zeros_like(turbine.af_data)

        for i,section in enumerate(turbine.af_data):
            # assume airfoil section as AOA of zero for slope calculations - for now
            a0_ind = section[0]['Alpha'].index(np.min(np.abs(section[0]['Alpha'])))
            # Coefficients
            if section[0]['NumTabs'] == 3:  # sections with flaps
                Clm[i,] = section[0]['Cl'][a0_ind]
                Cdm[i,] = section[0]['Cd'][a0_ind]
                Cl0[i,] = section[1]['Cl'][a0_ind]
                Cd0[i,] = section[1]['Cd'][a0_ind]
                Clp[i,] = section[2]['Cl'][a0_ind]
                Cdp[i,] = section[2]['Cd'][a0_ind]
                Ctrl_flp = float(section[2]['Ctrl'])
            else:                           # sections without flaps
                Cl0[i,] = Clp[i,] = Clm[i,] = section[0]['Cl'][a0_ind]
                Cd0[i,] = Cdp[i,] = Cdm[i,] = section[0]['Cd'][a0_ind]
                Ctrl = float(section[0]['Ctrl'])

        # Find slopes
        Kcl = (Clp - Cl0)/( (Ctrl_flp-Ctrl)*deg2rad )
        Kcd = (Cdp - Cd0)/( (Ctrl_flp-Ctrl)*deg2rad )

        # Find integrated constants
        kappa = np.zeros(len(v_rel))
        C1 = np.zeros(len(v_rel))
        C2 = np.zeros(len(v_rel))
        for i, (v_sec,phi) in enumerate(zip(v_rel, phi_vec)):
            C1[i] = integrate.trapz(0.5 * turbine.rho * turbine.chord * v_sec[0]**2 * turbine.span * Kcl * np.cos(phi))
            C2[i] = integrate.trapz(0.5 * turbine.rho * turbine.chord * v_sec[0]**2 * turbine.span * Kcd * np.sin(phi))
            kappa[i]=C1[i]+C2[i]

        # ------ Controller tuning -------
        # Open loop blade response
        zetaf  = turbine.bld_flapwise_damp
        omegaf = turbine.bld_flapwise_freq

        # Desired Closed loop response
        # zeta  = self.zeta_flp
        # omega = 4.6/(ts*zeta)

        # PI Gains
        if (self.zeta_flp == 0 or self.omega_flp == 0) or (not self.zeta_flp or not self.omega_flp):
            sys.exit('ERROR! --- Zeta and Omega flap must be nonzero for Flp_Mode >= 1 ---')

        self.Kp_flap = (2*self.zeta_flp*self.omega_flp - 2*zetaf*omegaf)/(kappa*omegaf**2)
        self.Ki_flap = (self.omega_flp**2 - omegaf**2)/(kappa*omegaf**2)

class ControllerBlocks():
    '''
    Class ControllerBlocks defines tuning parameters for additional controller features or "blocks"

    Methods:
    --------
    peak_shaving

    '''
    def __init__(self):
        pass

    def peak_shaving(self,controller, turbine):
        '''
        Define minimum blade pitch angle for peak shaving routine based on a maximum allowable thrust

        Parameters:
        -----------
        controller: class
                    Controller class containing controller operational information
        turbine: class
                 Turbine class containing necessary wind turbine information for controller tuning
        '''

        # Re-define Turbine Parameters for shorthand
        J = turbine.J                           # Total rotor inertial (kg-m^2)
        rho = turbine.rho                       # Air density (kg/m^3)
        R = turbine.rotor_radius                    # Rotor radius (m)
        A = np.pi*R**2                         # Rotor area (m^2)
        Ng = turbine.Ng                         # Gearbox ratio (-)
        rated_rotor_speed = turbine.rated_rotor_speed               # Rated rotor speed (rad/s)

        # Initialize some arrays
        Ct_op = np.empty(len(controller.TSR_op),dtype='float64')
        Ct_max = np.empty(len(controller.TSR_op),dtype='float64')
        beta_min = np.empty(len(controller.TSR_op),dtype='float64')
        # Find unshaved rotor thurst coefficients and associated rotor thrusts
        # for i in len(controller.TSR_op):
        for i in range(len(controller.TSR_op)):
            Ct_op[i] = turbine.Ct.interp_surface(controller.pitch_op[i],controller.TSR_op[i])
            T = 0.5 * rho * A * controller.v**2 * Ct_op

        # Define minimum max thrust and initialize pitch_min
        Tmax = controller.ps_percent * np.max(T)
        pitch_min = np.ones(len(controller.pitch_op)) * controller.min_pitch

        # Modify pitch_min if max thrust exceeds limits
        for i in range(len(controller.TSR_op)):
            # Find Ct values for operational TSR
            # Ct_tsr = turbine.Ct.interp_surface(turbine.pitch_initial_rad, controller.TSR_op[i])
            Ct_tsr = turbine.Ct.interp_surface(turbine.pitch_initial_rad,controller.TSR_op[i])
            # Define max Ct values
            Ct_max[i] = Tmax/(0.5 * rho * A * controller.v[i]**2)
            if T[i] > Tmax:
                Ct_op[i] = Ct_max[i]
            else:
                Ct_max[i] = np.minimum( np.max(Ct_tsr), Ct_max[i])
            # Define minimum pitch angle
            f_pitch_min = interpolate.interp1d(Ct_tsr, turbine.pitch_initial_rad, bounds_error=False, fill_value=(turbine.pitch_initial_rad[0],turbine.pitch_initial_rad[-1]))
            pitch_min[i] = f_pitch_min(Ct_max[i])

        controller.ps_min_bld_pitch = pitch_min

        # save some outputs for analysis or future work
        self.Tshaved = 0.5 * rho * A * controller.v**2 * Ct_op
        self.pitch_min = pitch_min
        self.v = controller.v
        self.Ct_max = Ct_max
        self.Ct_op = Ct_op
        self.T = T

    def min_pitch_saturation(self, controller, turbine):

        # Find TSR associated with minimum rotor speed
        TSR_at_minspeed = (controller.pc_minspd/turbine.Ng) * turbine.rotor_radius / controller.v_below_rated
        for i in range(len(TSR_at_minspeed)):
            if TSR_at_minspeed[i] > controller.TSR_op[i]:
                controller.TSR_op[i] = TSR_at_minspeed[i]

                # Initialize some arrays
                Cp_op = np.empty(len(turbine.pitch_initial_rad),dtype='float64')
                min_pitch = np.empty(len(TSR_at_minspeed),dtype='float64')

                # Find Cp-maximizing minimum pitch schedule
                # Find Cp coefficients at below-rated tip speed ratios
                Cp_op = turbine.Cp.interp_surface(turbine.pitch_initial_rad,TSR_at_minspeed[i])
                Cp_max = max(Cp_op)
                f_pitch_min = interpolate.interp1d(Cp_op, turbine.pitch_initial_rad, bounds_error=False, fill_value=(turbine.pitch_initial_rad[0],turbine.pitch_initial_rad[-1]))
                min_pitch[i] = f_pitch_min(Cp_max)

                # modify existing minimum pitch schedule
                controller.ps_min_bld_pitch[i] = np.maximum(controller.ps_min_bld_pitch[i], min_pitch[i])
            else:
                return

class ControllerTypes():
    '''
    Class ControllerTypes used to define any types of controllers that can be tuned.
        Generally, calculates gains based on some pre-defined tuning parameters.

    Methods:
    --------
    second_order_PI
    '''
    def __init__(self):
        pass

    def second_order_PI(self,zeta,om_n,A,B,linearize=False,v=None):
        '''
        Define proportional integral gain schedule for a closed
            loop system with a standard second-order form.

        Parameters:
        -----------
        zeta : int (-)
               Desired damping ratio
        om_n : int (rad/s)
               Desired natural frequency
        A : array_like (1/s)
            Plant poles (state transition matrix)
        B : array_like (varies)
            Plant numerators (input matrix)
        linearize : bool, optional
                    If 'True', find a gain scheduled based on a linearized plant.
        v : array_like (m/s)
            Wind speeds for linearized plant model, if desired.
        '''
        # Linearize system coefficients w.r.t. wind speed if desired
        if linearize:
            pA = np.polyfit(v,A,1)
            pB = np.polyfit(v,B,1)
            A = pA[0]*v + pA[1]
            B = pB[0]*v + pB[1]

        # Calculate gain schedule
        self.Kp = 1/B * (2*zeta*om_n + A)
        self.Ki = om_n**2/B


text:
import math
import unittest

import numpy as np
import numpy.testing as npt

from . import vector_utils as vu

class TestVectorUtils(unittest.TestCase):
    tol = 0.0001

    def test_normalize(self):
        vector = np.asarray([1, 1, 1, 1], dtype=np.float32)
        npt.assert_allclose([0.5, 0.5, 0.5, 0.5], vu.normalize(vector), rtol=self.tol)

    def test_cosine_similarity(self):
        npt.assert_almost_equal(vu.cosine_similarity([1, 0], [0, 1]), 0)
        npt.assert_almost_equal(vu.cosine_similarity([1, 1], [0, 1]), math.sqrt(2)/2)

    def test_cosine_similarity_complex(self):
        npt.assert_almost_equal(vu.cosine_similarity([1j], [1j]), 1)
        npt.assert_almost_equal(vu.cosine_similarity([1j], [1 + 1j]), math.sqrt(2)/2)
        npt.assert_almost_equal(vu.cosine_similarity([1 + 1j, 0], [1j, 1]), 0.5)
        npt.assert_almost_equal(vu.cosine_similarity([1, 1j], [1j, 1]), 0)

    def test_circular_convolution(self):
        npt.assert_allclose([1, 1], vu.circular_convolution([1, 0], [1, 1]),  rtol=self.tol)
        npt.assert_allclose([-1., 2., -1.], vu.circular_convolution([0, -1, 1], [1, 2, 3]),  rtol=self.tol)

    def test_create_vector(self):
        vec = vu.create_dense_random_vector(5, seed=2)
        npt.assert_allclose([-0.1280102, -0.94814754, 0.09932496, -0.12935521, -0.1592644], vec, rtol=self.tol)

        # Do it again to be sure.
        vec2 = vu.create_dense_random_vector(5, seed=2)
        npt.assert_allclose(vec, vec2, rtol=self.tol)

        # Make sure a different seed gives a different answer.
        vec3 = vu.create_dense_random_vector(5, seed=3)
        npt.assert_allclose([ 0.101596,  0.416296, -0.418191,  0.021655,  0.785894], vec3, rtol=self.tol)

    def test_create_complex_vector(self):
        vec = vu.create_dense_random_vector(1, seed=2, field=np.complex)
        self.assertTrue(np.issubdtype(vec.dtype, np.complex128))
        self.assertFalse(np.issubdtype(vec.dtype, np.float64))
        npt.assert_almost_equal(-0.1280102 - 0.94814754j, vec[0], self.tol)

    def test_get_k_neighbors_from_pairs(self):
        pairs = [('x', (1, 0, 0)), ('y', (0, 1, 0)), ('z', (0, 0, 1))]
        nearest = vu.get_k_neighbors_from_pairs(pairs, (0.9, 0.2, 0.1), 2)
        self.assertListEqual([nearest[0][0], nearest[1][0]], ['x', 'y'])
        npt.assert_allclose([0.971, 0.216], [nearest[0][1], nearest[1][1]], rtol=0.01)

        # Check no overflow
        nearest = vu.get_k_neighbors_from_pairs(pairs, (0.9, 0.2, 0.1), 10)
        self.assertIsNotNone(nearest)

    def test_complex_normalize(self):
        vector = np.asarray([1, 1j, -1, -1j], dtype=np.complex)
        npt.assert_allclose([0.5, 0.5j, -0.5, -0.5j], vu.normalize(vector), rtol=self.tol)

max_stars_repo_name: semanticvectors/semvecpy

max_issues_repo_name: semanticvectors/semvecpy

max_forks_repo_name: semanticvectors/semvecpy


text:
#!/usr/bin/env python3

import numpy as np
import os
import sys

ROOT_DIR = os.path.abspath("../")
sys.path.append(ROOT_DIR)
from utils import base

def laplacian(W):
    '''
    Construct the (weighted) Laplacian matrix from the (weighted) adjacnecy matrix

    Arguments:
    1. W:    weighted adjacency matrix

    Returns:
    1. L:    weighted Laplacian matrix
    '''
    assert type(W) == np.ndarray, "W must be of type 'numpy.ndarray'"
    assert W.size > 0, "W must not be empty"
    assert W.dtype == int or W.dtype == float, "W must be of dtype 'int' or 'float'"
    assert np.isfinite(W).all(), "Elements of W must be finite real numbers"
    size = W.shape
    assert len(size) == 2, "W must be 2D shape"
    assert size[0] == size[1], "W must be a square matrix"
    assert (np.diag(W) == 0).all(), "W must not have self-loop"

    # Construct the (weighted) Laplacian matrix
    # NOTE: Be-careful when W is a directed network
    #       in-link and out-link can be different
    L = np.diag(np.sum(W, 1)) - W

    return L

def hidden_effect(W, measure_id, hidden_id, a=0):
    '''
    Compute hidden node effect (C matrix)

    Arguments:
    1. W:               Weighted adjacency matrix
    2. measure_id:      Measured node indices
    3. hidden_id:       Hidden node indices
    4. a:               Dynamical constant -f'(X0)  (default: 0)

    Returns:
    1. C:       C matrix
    '''
    assert type(W) == np.ndarray, "W must be of type 'numpy.ndarray'"
    assert W.size > 0, "W must not be empty"
    assert W.dtype == int or W.dtype == float, "W must be of dtype 'int' or 'float'"
    assert np.isfinite(W).all(), "Elements of W must be finite real numbers"
    size = W.shape
    assert len(size) == 2, "W must be 2D shape"
    assert size[0] == size[1], "W must be a square matrix"
    assert (np.diag(W) == 0).all(), "W must not have self-loop"
    assert np.allclose(W, W.T), "W must be symmetric\n C matrix currently is defined for bi-directional network only"

    assert type(measure_id) == np.ndarray, "measure_id must be of type 'np.ndarray'"
    assert measure_id.size > 0, "measure_id must not be empty"
    assert measure_id.dtype == int, "measure_id must be of dtype 'int'"
    assert len(measure_id.shape) == 1, "measure_id must be 1D shape"
    assert (measure_id >= 0).all(), "measure_id elements must be non-negative integers"
    assert np.max(measure_id) < size[0], "measure_id elements must be smaller than the input matrix size"

    assert type(hidden_id) == np.ndarray, "hidden_id must be of type 'np.ndarray'"
    assert hidden_id.size > 0, "hidden_id must not be empty"
    assert hidden_id.dtype == int, "hidden_id must be of dtype 'int'"
    assert (hidden_id >= 0).all(), "hidden_id elements must be non-negative integers"
    assert len(hidden_id.shape) == 1, "hidden_id must be 1D shape"
    assert np.max(hidden_id) < size[0], "hidden_id elements must be smaller than the input matrix size"

    num_m = len(measure_id)
    num_h = len(hidden_id)
    assert size[0] == (num_m + num_h), "Total number of elements in measure_id and hidden_id does not equal the number of rows of M"

    all_id = np.unique(np.concatenate((measure_id, hidden_id)))
    assert len(all_id) == size[0], "mesure_id and hidden_id contain common elements"

    assert (type(a) == int or type(a) == float) and np.isfinite(a), "Dynamical constant a must be a real number"

    # Compute the Laplacian matrix
    L = laplacian(W)

    # Re-arrange the Laplacian matrix
    E = base.block_off_up(L, measure_id, hidden_id)
    Lh = base.block_diag_low(L, hidden_id)

    # Compute the C matrix
    H = Lh + a * np.identity(num_h)
    H_inv = base.inverse(H)
    C = np.matmul(np.matmul(E, H_inv), E.T)

    if not np.allclose(C, C.T):
        print("[WARN] The computed C matrix is not symmetric")

    return C

max_stars_repo_name: newTypeGeek/Network-Reconstruction

max_issues_repo_name: newTypeGeek/Network-Reconstruction

max_forks_repo_name: newTypeGeek/Network-Reconstruction


text:
import numpy as np

def centeredfft2(Z, FsX, FsY):
    """
    Computes a 2D fft centered at 0. Or in other words, compute
    the 2D fft, and then fftshift it.
    """
    M = np.size(Z, 1)
    N = np.size(Z, 0)

    k = np.array(range(-M // 2, M // 2))
    freqX = k / (M / FsX)

    k = np.array(range(-N // 2, N // 2))
    freqY = k / (N / FsY)

    fft = np.fft.fft2(Z) / (M * N)
    fft2 = (1 / (FsY * FsX)) * np.fft.fft2(Z)
    fft = np.fft.fftshift(fft)
    fft2 = np.fft.fftshift(fft2)

    return (freqX, freqY, fft, fft2)

max_stars_repo_name: collielimabean/LaserPolishingModel

max_issues_repo_name: collielimabean/LaserPolishingModel

max_forks_repo_name: collielimabean/LaserPolishingModel