text:
# Ref: http://francescopochetti.com/stock-market-prediction-part-ii-feature-generation/
from sklearn import preprocessing
import pandas_datareader.data as web
import pandas as pd
import classifier
import numpy as np
import operator

pd.options.mode.chained_assignment = None

def preprocessData(dataset):

    le = preprocessing.LabelEncoder()

    # in case divid-by-zero
    dataset.Open[dataset.Open == 0] = 1

    # add prediction target: next day Up/Down
    threshold = 0.000
    dataset['UpDown'] = (dataset['Close'] - dataset['Open']) / dataset['Open']
    dataset.UpDown[dataset.UpDown >= threshold] = 'Up'
    dataset.UpDown[dataset.UpDown < threshold] = 'Down'
    dataset.UpDown = le.fit(dataset.UpDown).transform(dataset.UpDown)
    dataset.UpDown = dataset.UpDown.shift(-1) # shift 1, so the y is actually next day's up/down
    dataset = dataset.drop(dataset.index[-1]) # drop last one because it has no up/down value
    return dataset

def count_missing(dataframe):
    return (dataframe.shape[0] * dataframe.shape[1]) - dataframe.count().sum()

def performCV(X_train, y_train, folds, method, parameters, savemodel):
    """
    given complete dataframe, number of folds, the % split to generate
    train and test set and features to perform prediction --> splits
    dataframein test and train set. Takes train set and splits in k folds.
    - Train on fold 1, test on 2
    - Train on fold 1-2, test on 3
    - Train on fold 1-2-3, test on 4
    ....
    returns mean of test accuracies
    """
    print ''
    print 'Parameters --------------------------------> ', parameters
    print 'Size train set: ', X_train.shape

    k = int(np.floor(float(X_train.shape[0])/folds))

    print 'Size of each fold: ', k

    acc = np.zeros(folds-1)
    for i in range(2, folds+1):
        print ''
        split = float(i-1)/i
        print 'Splitting the first ' + str(i) + ' chuncks at ' + str(i-1) + '/' + str(i)
        data = X_train[:(k*i)]
        output = y_train[:(k*i)]
        print 'Size of train+test: ', data.shape
        index = int(np.floor(data.shape[0]*split))
        X_tr = data[:index]
        y_tr = output[:index]

        X_te = data[(index+1):]
        y_te = output[(index+1):]

        acc[i-2] = classifier.performClassification(X_tr, y_tr, X_te, y_te, method, parameters, savemodel)
        print 'Accuracy on fold ' + str(i) + ': ', acc[i-2]

    return acc.mean()

def performTimeSeriesSearchGrid(X_train, y_train, folds, method, grid, savemodel):
    """
    parameters is a dictionary with: keys --> parameter , values --> list of values of parameter
    """
    print ''
    print 'Performing Search Grid CV...'
    print 'Algorithm: ', method
    param = grid.keys()
    finalGrid = {}
    if len(param) == 1:
        for value_0 in grid[param[0]]:
            parameters = [value_0]
            accuracy = performCV(X_train, y_train, folds, method, parameters, savemodel)
            finalGrid[accuracy] = parameters
        final = sorted(finalGrid.iteritems(), key=operator.itemgetter(0), reverse=True)
        print ''
        print finalGrid
        print ''
        print 'Final CV Results: ', final
        return final[0]

    elif len(param) == 2:
        for value_0 in grid[param[0]]:
            for value_1 in grid[param[1]]:
                parameters = [value_0, value_1]
                accuracy = performCV(X_train, y_train, folds, method, parameters, savemodel)
                finalGrid[accuracy] = parameters
        final = sorted(finalGrid.iteritems(), key=operator.itemgetter(0), reverse=True)
        print ''
        print finalGrid
        print ''
        print 'Final CV Results: ', final
        return final[0]

def performFeatureSelection(stock_name, maxdeltas, start, end, start_test, savemodel, method, folds, parameters):
    """
    """
    finalGrid = {}
    print ''
    print '============================================================='
    print ''
    for maxdelta in range(3, maxdeltas, 2):
        dataset = get_data(stock_name, start, end)
        delta = range(1, maxdelta)
        print 'Delta days accounted: ', max(delta)
        dataset = applyFeatures(dataset, delta)
        dataset = preprocessData(dataset)
        print 'Number of NaN: ', count_missing(dataset)
        X_train, y_train, X_test, y_test  = \
            classifier.prepareDataForClassification(dataset, start_test)

        print ''
        accuracy = performCV(X_train, y_train, folds, method, parameters, savemodel)
        finalGrid[accuracy] = maxdelta

    final = sorted(finalGrid.iteritems(), key=operator.itemgetter(0), reverse=True)
    print ''
    print finalGrid
    print ''
    print 'Final CV Results: ', final
    return final[0]

def performParameterSelection(stock_name, bestdelta, start, end, start_test, savemodel, method, folds, grid):
    """
    """
    dataset = get_data(stock_name, start, end)
    delta = range(1, bestdelta + 1)
    print 'Delta days accounted: ', max(delta)
    dataset = applyFeatures(dataset, delta)
    dataset = preprocessData(dataset)
    X_train, y_train, X_test, y_test  = \
            classifier.prepareDataForClassification(dataset, start_test)

    return performTimeSeriesSearchGrid(X_train, y_train, folds, method, grid, savemodel)

def addFeatures(dataframe, close, returns, n):
    """
    operates on two columns of dataframe:
    - append previous n days' OHLC and Volumn information
    - given Return_* computes the return of day i respect to day i-n.
    - given RolMean_* computes its moving average on n days

    """
    for c in dataframe.columns[0:5]:
        dataframe[c + str(n)] = dataframe[c].shift(n)

    return_n = "Return" + str(n)
    dataframe[return_n] = dataframe[close].pct_change(n)

    roll_n = "RolMean" + str(n)
    dataframe[roll_n] = dataframe[returns].rolling(n).mean()

def applyFeatures(dataset, delta):
    """
    applies rolling mean and delayed returns to each dataframe in the list
    """
    columns = dataset.columns
    close = columns[-3]
    returns = columns[-1]
    for n in delta:
        addFeatures(dataset, close, returns, n)

    dataset = dataset.drop(dataset.index[0:max(delta)]) #drop NaN due to delta spanning

    # normalize columns
    scaler = preprocessing.MinMaxScaler()
    return pd.DataFrame(scaler.fit_transform(dataset),\
            columns=dataset.columns, index=dataset.index)

def get_data(name, start, end):
    data = web.get_data_yahoo(name, start, end)
    del data['Adj Close'] # we don't need Adj Close
    data['Return'] = (data['Close'] - data['Open']) / data['Open']

    return data

max_stars_repo_name: rbcwmlab/stock-price-prediction

max_issues_repo_name: Rgveda/stock-price-prediction

max_forks_repo_name: Rgveda/stock-price-prediction


text:
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

train_data = pd.read_csv("Chapter01/mtcars.csv")

X = np.array(train_data.hp).reshape(32,1)
y = np.array(train_data.mpg).reshape(32,1)

model = LinearRegression(fit_intercept=True)
model.fit(X, y)

print(model.intercept_[0], model.coef_[0])

def gradient_descent(x, y, learning_rate, conv_threshold, batch_size, max_iter):
    converged = False
    iteration = 0
    m = batch_size
    t0 = np.random.random(x.shape[1])
    t1 = np.random.random(x.shape[1])
    MSE = (sum([(t0 + t1*x[i] - y[i])**2 for i in range(m)])) / m
    while not converged:
        grad0 = 1.0 / m * sum([(t0 + t1*x[i] - y[i]) for i in range(m)]) # theta0 미분
        grad1 = 1.0 / m * sum([(t0 + t1*x[i] - y[i]) * x[i] for i in range(m)]) # theta1 미분
        temp0 = t0 - learning_rate * grad0
        temp1 = t1 - learning_rate * grad1
        t0 = temp0
        t1 = temp1
        MSE_new = (sum([(t0 + t1*x[i] - y[i])**2 for i in range(m)])) / m
        if iteration % 200 == 0:
            print(MSE_new)
        if abs(MSE - MSE_new) <= conv_threshold:
            print(iteration)
            converged = True
        MSE = MSE_new
        iteration += 1
        if iteration == max_iter:
            converged = True
    return t0, t1

# main
Inter, Coeff = gradient_descent(x = X,y = y, learning_rate=0.00003,conv_threshold=1e-8, batch_size=32,max_iter=1500000)
print(Inter, Coeff)

ext: py

max_stars_repo_name: cutz-j/Statistics-for-ML

max_issues_repo_name: cutz-j/Statistics-for-ML

max_forks_repo_name: cutz-j/Statistics-for-ML

reason: import numpy


text:
# Copyright 2018 Amazon.com, Inc. or its affiliates. All Rights Reserved.
#
#   Licensed under the Apache License, Version 2.0 (the "License").
#   You may not use this file except in compliance with the License.
#   A copy of the License is located at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   or in the "license" file accompanying this file. This file 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 pytest
import mxnet as mx
import numpy as np
from mxfusion.components.variables.runtime_variable import add_sample_dimension, array_has_samples, get_num_samples
from mxfusion.components.distributions import Normal, MultivariateNormal
from mxfusion.util.testutils import numpy_array_reshape, plot_univariate, plot_bivariate
from mxfusion.util.testutils import MockMXNetRandomGenerator

@pytest.mark.usefixtures("set_seed")
class TestNormalDistribution(object):

    @pytest.mark.parametrize("dtype, mean, mean_isSamples, var, var_isSamples, rv, rv_isSamples, num_samples", [
        (np.float64, np.random.rand(5,3,2), True, np.random.rand(3,2)+0.1, False, np.random.rand(5,3,2), True, 5),
        (np.float64, np.random.rand(3,2), False, np.random.rand(5,3,2)+0.1, True, np.random.rand(5,3,2), True, 5),
        (np.float64, np.random.rand(3,2), False, np.random.rand(3,2)+0.1, False, np.random.rand(5,3,2), True, 5),
        (np.float64, np.random.rand(3,2), False, np.random.rand(3,2)+0.1, False, np.random.rand(3,2), False, 1),
        (np.float32, np.random.rand(5,3,2), True, np.random.rand(3,2)+0.1, False, np.random.rand(5,3,2), True, 5),
        ])
    def test_log_pdf(self, dtype, mean, mean_isSamples, var, var_isSamples,
                     rv, rv_isSamples, num_samples):
        from scipy.stats import norm

        isSamples_any = any([mean_isSamples, var_isSamples, rv_isSamples])
        rv_shape = rv.shape[1:] if rv_isSamples else rv.shape
        n_dim = 1 + len(rv.shape) if isSamples_any and not rv_isSamples else len(rv.shape)
        mean_np = numpy_array_reshape(mean, mean_isSamples, n_dim)
        var_np = numpy_array_reshape(var, var_isSamples, n_dim)
        rv_np = numpy_array_reshape(rv, rv_isSamples, n_dim)
        log_pdf_np = norm.logpdf(rv_np, mean_np, np.sqrt(var_np))
        normal = Normal.define_variable(shape=rv_shape, dtype=dtype).factor
        mean_mx = mx.nd.array(mean, dtype=dtype)
        if not mean_isSamples:
            mean_mx = add_sample_dimension(mx.nd, mean_mx)
        var_mx = mx.nd.array(var, dtype=dtype)
        if not var_isSamples:
            var_mx = add_sample_dimension(mx.nd, var_mx)
        rv_mx = mx.nd.array(rv, dtype=dtype)
        if not rv_isSamples:
            rv_mx = add_sample_dimension(mx.nd, rv_mx)
        variables = {normal.mean.uuid: mean_mx, normal.variance.uuid: var_mx, normal.random_variable.uuid: rv_mx}
        log_pdf_rt = normal.log_pdf(F=mx.nd, variables=variables)

        assert np.issubdtype(log_pdf_rt.dtype, dtype)
        assert array_has_samples(mx.nd, log_pdf_rt) == isSamples_any
        if isSamples_any:
            assert get_num_samples(mx.nd, log_pdf_rt) == num_samples
        if np.issubdtype(dtype, np.float64):
            rtol, atol = 1e-7, 1e-10
        else:
            rtol, atol = 1e-4, 1e-5
        assert np.allclose(log_pdf_np, log_pdf_rt.asnumpy(), rtol=rtol, atol=atol)

    @pytest.mark.parametrize(
        "dtype, mean, mean_isSamples, var, var_isSamples, rv_shape, num_samples",[
        (np.float64, np.random.rand(5,3,2), True, np.random.rand(3,2)+0.1, False, (3,2), 5),
        (np.float64, np.random.rand(3,2), False, np.random.rand(5,3,2)+0.1, True, (3,2), 5),
        (np.float64, np.random.rand(3,2), False, np.random.rand(3,2)+0.1, False, (3,2), 5),
        (np.float64, np.random.rand(5,3,2), True, np.random.rand(5,3,2)+0.1, True, (3,2), 5),
        (np.float32, np.random.rand(5,3,2), True, np.random.rand(3,2)+0.1, False, (3,2), 5),
        ])
    def test_draw_samples(self, dtype, mean, mean_isSamples, var,
                          var_isSamples, rv_shape, num_samples):
        n_dim = 1 + len(rv_shape)
        mean_np = numpy_array_reshape(mean, mean_isSamples, n_dim)
        var_np = numpy_array_reshape(var, var_isSamples, n_dim)

        rand = np.random.randn(num_samples, *rv_shape)
        rv_samples_np = mean_np + rand * np.sqrt(var_np)

        rand_gen = MockMXNetRandomGenerator(mx.nd.array(rand.flatten(), dtype=dtype))

        normal = Normal.define_variable(shape=rv_shape, dtype=dtype,
                                        rand_gen=rand_gen).factor
        mean_mx = mx.nd.array(mean, dtype=dtype)
        if not mean_isSamples:
            mean_mx = add_sample_dimension(mx.nd, mean_mx)
        var_mx = mx.nd.array(var, dtype=dtype)
        if not var_isSamples:
            var_mx = add_sample_dimension(mx.nd, var_mx)
        variables = {normal.mean.uuid: mean_mx, normal.variance.uuid: var_mx}
        rv_samples_rt = normal.draw_samples(
            F=mx.nd, variables=variables, num_samples=num_samples)

        assert np.issubdtype(rv_samples_rt.dtype, dtype)
        assert array_has_samples(mx.nd, rv_samples_rt)
        assert get_num_samples(mx.nd, rv_samples_rt) == num_samples

        if np.issubdtype(dtype, np.float64):
            rtol, atol = 1e-7, 1e-10
        else:
            rtol, atol = 1e-4, 1e-5
        assert np.allclose(rv_samples_np, rv_samples_rt.asnumpy(), rtol=rtol, atol=atol)

    def test_draw_samples_non_mock(self, plot=False):
        # Also make sure the non-mock sampler works
        dtype = np.float32
        num_samples = 100000

        mean = np.array([0.5])
        variance = np.array([2])

        rv_shape = (1,)

        mean_mx = add_sample_dimension(mx.nd, mx.nd.array(mean, dtype=dtype))
        variance_mx = add_sample_dimension(mx.nd, mx.nd.array(variance, dtype=dtype))

        rand_gen = None
        var = Normal.define_variable(shape=rv_shape, rand_gen=rand_gen, dtype=dtype).factor
        variables = {var.mean.uuid: mean_mx, var.variance.uuid: variance_mx}
        rv_samples_rt = var.draw_samples(F=mx.nd, variables=variables, num_samples=num_samples)

        assert array_has_samples(mx.nd, rv_samples_rt)
        assert get_num_samples(mx.nd, rv_samples_rt) == num_samples
        assert rv_samples_rt.dtype == dtype

        from scipy.stats import norm
        if plot:
            plot_univariate(samples=rv_samples_rt, dist=norm, loc=mean[0], scale=np.sqrt(variance[0]))

        mean_est, scale_est = norm.fit(rv_samples_rt.asnumpy().ravel())
        mean_tol = 1e-1
        variance_tol = 1e-1
        assert np.abs(mean[0] - mean_est) < mean_tol
        assert np.abs(variance[0] - scale_est ** 2) < variance_tol

def make_symmetric(array):
    original_shape = array.shape
    d3_array = np.reshape(array, (-1,)+array.shape[-2:])
    d3_array = (d3_array[:,:,:,None]*d3_array[:,:,None,:]).sum(-3)+np.eye(2)
    return np.reshape(d3_array, original_shape)

@pytest.mark.usefixtures("set_seed")
class TestMultivariateNormalDistribution(object):

    @pytest.mark.parametrize("dtype, mean, mean_isSamples, var, var_isSamples, rv, rv_isSamples, num_samples", [
        (np.float32, np.random.rand(3,2), False, make_symmetric(np.random.rand(5,3,2,2)+0.1), True, np.random.rand(5,3,2), True, 5),
        (np.float32, np.random.rand(5,3,2), True, make_symmetric(np.random.rand(3,2,2)+0.1), False, np.random.rand(5,3,2), True, 5),
        (np.float32, np.random.rand(3,2), False, make_symmetric(np.random.rand(3,2,2)+0.1), False, np.random.rand(5,3,2), True, 5),
        ])
    def test_log_pdf_with_broadcast(self, dtype, mean, mean_isSamples, var, var_isSamples,
                        rv, rv_isSamples, num_samples):

        mean_mx = mx.nd.array(mean, dtype=dtype)
        if not mean_isSamples:
            mean_mx = add_sample_dimension(mx.nd, mean_mx)
        mean = mean_mx.asnumpy()

        var_mx = mx.nd.array(var, dtype=dtype)
        if not var_isSamples:
            var_mx = add_sample_dimension(mx.nd, var_mx)
        var = var_mx.asnumpy()

        rv_mx = mx.nd.array(rv, dtype=dtype)
        if not rv_isSamples:
            rv_mx = add_sample_dimension(mx.nd, rv_mx)
        rv = rv_mx.asnumpy()

        print(mean_mx.shape, var_mx.shape, rv_mx.shape)

        from scipy.stats import multivariate_normal
        isSamples_any = any([mean_isSamples, var_isSamples, rv_isSamples])
        rv_shape = rv.shape[1:]

        n_dim = 1 + len(rv.shape) if isSamples_any and not rv_isSamples else len(rv.shape)
        mean_np = np.broadcast_to(mean,(5,3,2))
        var_np = np.broadcast_to(var,(5,3,2,2))
        rv_np = numpy_array_reshape(rv, isSamples_any, n_dim)

        rand = np.random.rand(num_samples, *rv_shape)
        rand_gen = MockMXNetRandomGenerator(mx.nd.array(rand.flatten(), dtype=dtype))

        r = []
        for s in range(len(rv_np)):
            a = []
            for i in range(len(rv_np[s])):
                a.append(multivariate_normal.logpdf(rv_np[s][i], mean_np[s][i], var_np[s][i]))
            r.append(a)
        log_pdf_np = np.array(r)

        normal = MultivariateNormal.define_variable(shape=rv_shape, dtype=dtype, rand_gen=rand_gen).factor
        variables = {normal.mean.uuid: mean_mx, normal.covariance.uuid: var_mx, normal.random_variable.uuid: rv_mx}
        log_pdf_rt = normal.log_pdf(F=mx.nd, variables=variables)

        assert np.issubdtype(log_pdf_rt.dtype, dtype)
        assert array_has_samples(mx.nd, log_pdf_rt) == isSamples_any
        if isSamples_any:
            assert get_num_samples(mx.nd, log_pdf_rt) == num_samples, (get_num_samples(mx.nd, log_pdf_rt), num_samples)
        assert np.allclose(log_pdf_np, log_pdf_rt.asnumpy())

    @pytest.mark.parametrize("dtype, mean, mean_isSamples, var, var_isSamples, rv, rv_isSamples, num_samples", [
        (np.float64, np.random.rand(5,3,2), True, make_symmetric(np.random.rand(5,3,2,2)+0.1), True, np.random.rand(5,3,2), True, 5),
        ])
    def test_log_pdf_no_broadcast(self, dtype, mean, mean_isSamples, var, var_isSamples,
                        rv, rv_isSamples, num_samples):

        mean_mx = mx.nd.array(mean, dtype=dtype)
        if not mean_isSamples:
            mean_mx = add_sample_dimension(mx.nd, mean_mx)
        mean = mean_mx.asnumpy()

        var_mx = mx.nd.array(var, dtype=dtype)
        if not var_isSamples:
            var_mx = add_sample_dimension(mx.nd, var_mx)
        var = var_mx.asnumpy()

        rv_mx = mx.nd.array(rv, dtype=dtype)
        if not rv_isSamples:
            rv_mx = add_sample_dimension(mx.nd, rv_mx)
        rv = rv_mx.asnumpy()

        from scipy.stats import multivariate_normal
        isSamples_any = any([mean_isSamples, var_isSamples, rv_isSamples])
        rv_shape = rv.shape[1:]

        n_dim = 1 + len(rv.shape) if isSamples_any and not rv_isSamples else len(rv.shape)
        mean_np = numpy_array_reshape(mean, isSamples_any, n_dim)
        var_np = numpy_array_reshape(var, isSamples_any, n_dim)
        rv_np = numpy_array_reshape(rv, isSamples_any, n_dim)

        rand = np.random.rand(num_samples, *rv_shape)
        rand_gen = MockMXNetRandomGenerator(mx.nd.array(rand.flatten(), dtype=dtype))

        r = []
        for s in range(len(rv_np)):
            a = []
            for i in range(len(rv_np[s])):
                a.append(multivariate_normal.logpdf(rv_np[s][i], mean_np[s][i], var_np[s][i]))
            r.append(a)
        log_pdf_np = np.array(r)

        normal = MultivariateNormal.define_variable(shape=rv_shape, dtype=dtype, rand_gen=rand_gen).factor
        variables = {normal.mean.uuid: mean_mx, normal.covariance.uuid: var_mx, normal.random_variable.uuid: rv_mx}
        log_pdf_rt = normal.log_pdf(F=mx.nd, variables=variables)

        assert np.issubdtype(log_pdf_rt.dtype, dtype)
        assert array_has_samples(mx.nd, log_pdf_rt) == isSamples_any
        if isSamples_any:
            assert get_num_samples(mx.nd, log_pdf_rt) == num_samples, (get_num_samples(mx.nd, log_pdf_rt), num_samples)
        assert np.allclose(log_pdf_np, log_pdf_rt.asnumpy())

    @pytest.mark.parametrize(
        "dtype, mean, mean_isSamples, var, var_isSamples, rv_shape, num_samples",[
        (np.float64, np.random.rand(3,2), False, make_symmetric(np.random.rand(3,2,2)+0.1), False, (3,2), 5),
        ])
    def test_draw_samples_with_broadcast(self, dtype, mean, mean_isSamples, var,
                             var_isSamples, rv_shape, num_samples):

        mean_mx = mx.nd.array(mean, dtype=dtype)
        if not mean_isSamples:
            mean_mx = add_sample_dimension(mx.nd, mean_mx)
        var_mx = mx.nd.array(var, dtype=dtype)
        if not var_isSamples:
            var_mx = add_sample_dimension(mx.nd, var_mx)
        var = var_mx.asnumpy()

        rand = np.random.rand(num_samples, *rv_shape)
        rand_gen = MockMXNetRandomGenerator(mx.nd.array(rand.flatten(), dtype=dtype))
        rv_samples_np = mean + np.matmul(np.linalg.cholesky(var), np.expand_dims(rand, axis=-1)).sum(-1)

        normal = MultivariateNormal.define_variable(shape=rv_shape, dtype=dtype, rand_gen=rand_gen).factor
        variables = {normal.mean.uuid: mean_mx, normal.covariance.uuid: var_mx}
        draw_samples_rt = normal.draw_samples(F=mx.nd, variables=variables, num_samples=num_samples)

        assert np.issubdtype(draw_samples_rt.dtype, dtype)
        assert np.allclose(rv_samples_np, draw_samples_rt.asnumpy())

    @pytest.mark.parametrize(
        "dtype, mean, mean_isSamples, var, var_isSamples, rv_shape, num_samples",[
        (np.float64, np.random.rand(3,2), False, make_symmetric(np.random.rand(3,2,2)+0.1), False, (3,2), 5),
        (np.float64, np.random.rand(3,2), False, make_symmetric(np.random.rand(5,3,2,2)+0.1), True, (3,2), 5),
        (np.float64, np.random.rand(5,3,2), True, make_symmetric(np.random.rand(3,2,2)+0.1), False, (3,2), 5),
        (np.float64, np.random.rand(5,3,2), True, make_symmetric(np.random.rand(5,3,2,2)+0.1), True, (3,2), 5),
        ])
    def test_draw_samples_no_broadcast(self, dtype, mean, mean_isSamples, var,
                             var_isSamples, rv_shape, num_samples):

        mean_mx = mx.nd.array(mean, dtype=dtype)
        if not mean_isSamples:
            mean_mx = add_sample_dimension(mx.nd, mean_mx)
        var_mx = mx.nd.array(var, dtype=dtype)
        if not var_isSamples:
            var_mx = add_sample_dimension(mx.nd, var_mx)
        var = var_mx.asnumpy()

        # n_dim = 1 + len(rv.shape) if isSamples_any else len(rv.shape)
        rand = np.random.rand(num_samples, *rv_shape)
        rand_gen = MockMXNetRandomGenerator(mx.nd.array(rand.flatten(), dtype=dtype))
        rand_exp = np.expand_dims(rand, axis=-1)
        lmat = np.linalg.cholesky(var)
        temp1 = np.matmul(lmat, rand_exp).sum(-1)
        rv_samples_np = mean + temp1

        normal = MultivariateNormal.define_variable(shape=rv_shape, dtype=dtype, rand_gen=rand_gen).factor

        variables = {normal.mean.uuid: mean_mx, normal.covariance.uuid: var_mx}
        draw_samples_rt = normal.draw_samples(F=mx.nd, variables=variables, num_samples=num_samples)

        assert np.issubdtype(draw_samples_rt.dtype, dtype)
        assert array_has_samples(mx.nd, draw_samples_rt) == True
        assert get_num_samples(mx.nd, draw_samples_rt) == num_samples, (get_num_samples(mx.nd, draw_samples_rt), num_samples)

    def test_draw_samples_non_mock(self, plot=False):
        # Also make sure the non-mock sampler works
        dtype = np.float32
        num_samples = 100000

        mean = np.array([0.5, 0])
        covariance = np.array([[2, 0.5], [0.5, 2]])

        rv_shape = (2,)

        mean_mx = add_sample_dimension(mx.nd, mx.nd.array(mean, dtype=dtype))
        covariance_mx = add_sample_dimension(mx.nd, mx.nd.array(covariance, dtype=dtype))

        rand_gen = None
        var = MultivariateNormal.define_variable(shape=rv_shape, rand_gen=rand_gen, dtype=dtype).factor
        variables = {var.mean.uuid: mean_mx, var.covariance.uuid: covariance_mx}
        rv_samples_rt = var.draw_samples(F=mx.nd, variables=variables, num_samples=num_samples)

        assert array_has_samples(mx.nd, rv_samples_rt)
        assert get_num_samples(mx.nd, rv_samples_rt) == num_samples
        assert rv_samples_rt.dtype == dtype

        from scipy.stats import multivariate_normal
        if plot:
            plot_bivariate(samples=rv_samples_rt, dist=multivariate_normal, mean=mean, cov=covariance)

        # mean_est, scale_est = multivariate_normal.fit(rv_samples_rt.asnumpy().ravel())
        mean_est = np.mean(rv_samples_rt.asnumpy(), axis=0)
        cov_est = np.cov(rv_samples_rt.asnumpy(), rowvar=False)
        mean_tol = 1e-1
        covariance_tol = 1e-1
        assert np.allclose(mean, mean_est, atol=mean_tol)
        assert np.allclose(covariance, cov_est, covariance_tol)

max_stars_repo_name: JeremiasKnoblauch/MXFusion

max_issues_repo_name: JeremiasKnoblauch/MXFusion

max_forks_repo_name: JeremiasKnoblauch/MXFusion

reason: import numpy,from scipy


text:
"""
Test script for testing the AiryBeam2D command.
"""

from LightPipes import *
import matplotlib.pyplot as plt
import numpy as np

wavelength = 2.3*um
size = 30*mm
N = 500
x0=y0=1*mm
a1=a2=0.1/mm
z=900*cm

F0=Begin(size,wavelength,N)
F0=AiryBeam2D(F0,x0=x0, y0=y0, a1=a1, a2=a2)
F=Fresnel(F0,z)
I=Intensity(F)
plt.figure(figsize = (9,5))
plt.imshow(I,
           extent=[-size/2/mm, size/2/mm, -size/2/mm, size/2/mm],
           origin='lower',
           cmap='jet',
           )
plt.title('2D Airy beam')
plt.xlabel('x [mm]')
plt.ylabel('y [mm]')
s = r'LightPipes for Python' + '\n'+ '2D Airy beam' + '\n\n'\
    r'$\lambda = {:4.2f}$'.format(wavelength/um) + r' ${\mu}m$' + '\n'\
    r'$size = {:4.2f}$'.format(size/mm) + r' $mm$' + '\n'\
    r'$N = {:4d}$'.format(N) + '\n'\
    r'$x_0 = y_0 = {:4.2f}$'.format(x0/mm) + r' $mm$' + '\n'\
    r'$a1 = a2 =  $' + '{:4.2f}'.format(a1*mm) + r' $/mm$' + '\n'\
    r'$z = $' + '{:4.2f}'.format(z/cm) + r' $cm$' + '\n'\
    r'${\copyright}$ Fred van Goor, May 2022'
plt.text(16, -10, s, bbox={'facecolor': 'white', 'pad': 5})
plt.show()

max_stars_repo_name: opticspy/lightpipes-python

max_issues_repo_name: opticspy/lightpipes-python

max_forks_repo_name: opticspy/lightpipes-python


text:
import numpy as np
from slfpy import SLF

def to_slf(obj,data):
  x=data['x']['data']
  y=data['y']['data']
  elem=data['elem']['data']
  dt=data['time']['data']
  slf=SLF()
  slf.addTITLE("")
  slf.addPOIN(np.column_stack((x,y)))
  slf.addIKLE(elem)
  minDate=np.min(dt)
  # TODO minDate
  seconds=(dt-minDate).astype("timedelta64[s]")
  slf.tags['times']=seconds

  del data['x']
  del data['y']
  del data['elem']
  del data['time']

  array=None
  for vname in data:
    slf.addVAR({"name":data[vname]['meta']['standard_name'],"unit":data[vname]['meta']['units']})
    if array is None:array=data[vname]['data']
    else:array=np.stack((array,data[vname]['data']))
  slf.values=np.einsum('ijk->jik', array)
  filepath=obj['filepath']+".slf"
  slf.write(filepath)
  return filepath


text:
import argparse
import os
import sys
import random
import timeit
import datetime
import numpy as np
import pickle
import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
import torch.backends.cudnn as cudnn
from torch.utils import data, model_zoo
from torch.autograd import Variable
import torchvision.transforms as transform
from model.deeplabv2 import Res_Deeplab
from utils.loss import CrossEntropy2d
from utils.loss import CrossEntropyLoss2dPixelWiseWeighted
from utils.loss import MSELoss2d
from utils import transformmasks
from utils import transformsgpu
from utils.helpers import colorize_mask
import utils.palette as palette
from utils.sync_batchnorm import convert_model
from utils.sync_batchnorm import DataParallelWithCallback
from data.voc_dataset import VOCDataSet
from data import get_loader, get_data_path
from data.augmentations import *
from tqdm import tqdm
import PIL
from torchvision import transforms
import json
from torch.utils import tensorboard
from evaluateSSL import evaluate
import time
from dl_lib.modeling.meta_arch import DynamicNet4Seg
from dl_lib.solver import build_optimizer
from dl_lib.checkpoint import DetectionCheckpointer
from dl_lib.engine import launch
from dl_lib.checkpoint import DetectionCheckpointer
import matplotlib.pyplot as plt
import matplotlib.image as mpimg  # mpimg 用于读取图片
import numpy as np

start = timeit.default_timer()
start_writeable = datetime.datetime.now().strftime('%m-%d_%H-%M')

def get_arguments():
    """Parse all the arguments provided from the CLI.

    Returns:
      A list of parsed arguments.
    """
    parser = argparse.ArgumentParser(description="DeepLab-ResNet Network")
    parser.add_argument("--gpus", type=int, default=1,
                        help="choose number of gpu devices to use (default: 1)")
    parser.add_argument("-c", "--config", type=str, default='config.json',
                        help='Path to the config file (default: config.json)')
    parser.add_argument("-r", "--resume", type=str, default=None,
                        help='Path to the .pth file to resume from (default: None)')
    parser.add_argument("-n", "--name", type=str, default=None, required=True,
                        help='Name of the run (default: None)')
    parser.add_argument("--save-images", type=str, default=None,
                        help='Include to save images (default: None)')
    return parser.parse_args()

def loss_calc(pred, label):
    label = Variable(label.long()).cuda()
    if len(gpus) > 1:
        criterion = torch.nn.DataParallel(CrossEntropy2d(ignore_label=ignore_label),
                                          device_ids=gpus).cuda()  # Ignore label ??
    else:
        criterion = CrossEntropy2d(ignore_label=ignore_label).cuda()  # Ignore label ??

    return criterion(pred, label)

def lr_poly(base_lr, iter, max_iter, power):
    return base_lr * ((1 - float(iter) / max_iter) ** (power))

def adjust_learning_rate(optimizer, i_iter):
    lr = lr_poly(learning_rate, i_iter, num_iterations, lr_power)
    optimizer.param_groups[0]['lr'] = lr
    if len(optimizer.param_groups) > 1:
        optimizer.param_groups[1]['lr'] = lr * 10

def sigmoid_ramp_up(iter, max_iter):
    if iter >= max_iter:
        return 1
    else:
        return np.exp(- 5 * (1 - iter / max_iter) ** 2)

def strongTransform(parameters, data=None, target=None):
    assert ((data is not None) or (target is not None))
    data, target = transformsgpu.mix(mask=parameters["Mix"], data=data, target=target)
    data, target = transformsgpu.colorJitter(colorJitter=parameters["ColorJitter"],
                                             img_mean=torch.from_numpy(IMG_MEAN.copy()).cuda(), data=data,
                                             target=target)
    data, target = transformsgpu.gaussian_blur(blur=parameters["GaussianBlur"], data=data, target=None)
    data, target = transformsgpu.flip(flip=parameters["flip"], data=data, target=target)
    return data, target

def weakTransform(parameters, data=None, target=None):
    data, target = transformsgpu.flip(flip=parameters["flip"], data=data, target=target)
    return data, target

def getWeakInverseTransformParameters(parameters):
    return parameters

def getStrongInverseTransformParameters(parameters):
    return parameters

class DeNormalize(object):
    def __init__(self, mean):
        self.mean = mean

    def __call__(self, tensor):
        IMG_MEAN = torch.from_numpy(self.mean.copy())
        IMG_MEAN, _ = torch.broadcast_tensors(IMG_MEAN.unsqueeze(1).unsqueeze(2), tensor)
        tensor = tensor + IMG_MEAN
        tensor = (tensor / 255).float()
        tensor = torch.flip(tensor, (0,))
        return tensor

class Learning_Rate_Object(object):
    def __init__(self, learning_rate):
        self.learning_rate = learning_rate

def save_image(image, epoch, id, palette):
    with torch.no_grad():
        if image.shape[0] == 3:
            restore_transform = transforms.Compose([
                DeNormalize(IMG_MEAN),
                transforms.ToPILImage()])

            image = restore_transform(image)
            # image = PIL.Image.fromarray(np.array(image)[:, :, ::-1])  # BGR->RGB
            image.save(os.path.join('../visualiseImages/', str(epoch) + id + '.png'))
        else:
            mask = image.numpy()
            colorized_mask = colorize_mask(mask, palette)
            colorized_mask.save(os.path.join('../visualiseImages/', str(epoch) + id + '.png'))

def _save_checkpoint(iteration, model, optimizer, config, ema_model, save_best=False, overwrite=True):
    checkpoint = {
        'iteration': iteration,
        'optimizer': optimizer.state_dict(),
        'config': config,
    }
    if len(gpus) > 1:
        checkpoint['model'] = model.module.state_dict()
    else:
        checkpoint['model'] = model.state_dict()

    if save_best:
        filename = os.path.join(checkpoint_dir, f'best_model.pth')
        torch.save(checkpoint, filename)
        print("Saving current best model: best_model.pth")
    else:
        filename = os.path.join(checkpoint_dir, f'checkpoint-iter{iteration}.pth')
        print(f'\nSaving a checkpoint: {filename} ...')
        torch.save(checkpoint, filename)
        if overwrite:
            try:
                os.remove(os.path.join(checkpoint_dir, f'checkpoint-iter{iteration - save_checkpoint_every}.pth'))
            except:
                pass

def _resume_checkpoint(resume_path, model, optimizer, ema_model):
    print(f'Loading checkpoint : {resume_path}')
    checkpoint = torch.load(resume_path)

    # Load last run info, the model params, the optimizer and the loggers
    iteration = checkpoint['iteration'] + 1
    print('Starting at iteration: ' + str(iteration))

    if len(gpus) > 1:
        model.module.load_state_dict(checkpoint['model'])
    else:
        model.load_state_dict(checkpoint['model'])

    optimizer.load_state_dict(checkpoint['optimizer'])

    return iteration, model, optimizer, ema_model

def main(config):
    print("batch_size:",config['training']['batch_size'],"lr:",config['SOLVER']['OPTIMIZER']['BASE_LR'])
    print(config)
    best_mIoU = 0

    if consistency_loss == 'CE':
        if len(gpus) > 1:
            unlabeled_loss = torch.nn.DataParallel(CrossEntropyLoss2dPixelWiseWeighted(ignore_index=ignore_label),
                                                   device_ids=gpus).cuda()
        else:
            unlabeled_loss = CrossEntropyLoss2dPixelWiseWeighted().cuda()
    elif consistency_loss == 'MSE':
        if len(gpus) > 1:
            unlabeled_loss = torch.nn.DataParallel(MSELoss2d(), device_ids=gpus).cuda()
        else:
            unlabeled_loss = MSELoss2d().cuda()
    cudnn.enabled = True

    model = DynamicNet4Seg(config)

    if len(gpus) > 1:  # parallel
        if use_sync_batchnorm:
            # port = 2 ** 15 + 2 ** 14 + hash(os.getuid()) % 2 ** 14
            import socket
            sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
            # Binding to port 0 will cause the OS to find an available port for us
            sock.bind(("", 0))
            port = sock.getsockname()[1]
            sock.close()
            dist_url = "tcp://127.0.0.1:{}".format(port)
            torch.distributed.init_process_group(backend="NCCL",
                                    init_method=dist_url,
                                    world_size=1,rank=0)
            # torch.distributed.init_process_group('nccl', init_method='file:///home/.../my_file', world_size=1, rank=0)
            model = convert_model(model)
            model = DataParallelWithCallback(model, device_ids=gpus)
        else:
            model = torch.nn.DataParallel(model, device_ids=gpus)
    model.train()
    model.cuda()
    cudnn.benchmark = True

    # dataset
    if dataset == 'pascal_voc':
        data_loader = get_loader(dataset)
        data_path = get_data_path(dataset)
        train_dataset = data_loader(data_path, crop_size=input_size, scale=random_scale, mirror=random_flip)
    elif dataset == 'cityscapes':
        data_loader = get_loader('cityscapes')
        data_path = get_data_path('cityscapes')
        if random_crop:
            data_aug = Compose([RandomCrop_city(input_size)])
        else:
            data_aug = None
        train_dataset = data_loader(data_path, is_transform=True, augmentations=data_aug, img_size=input_size)
    train_dataset_size = len(train_dataset)
    # unlabel dataset
    partial_size = labeled_samples
    print('Training on number of samples:', partial_size)
    if split_id is not None:
        train_ids = pickle.load(open(split_id, 'rb'))
        print('loading train ids from {}'.format(split_id))
    else:
        np.random.seed(random_seed)
        train_ids = np.arange(train_dataset_size)
        np.random.shuffle(train_ids)
    if train_unlabeled:
        train_sampler = data.sampler.SubsetRandomSampler(train_ids[:partial_size])
        trainloader = data.DataLoader(train_dataset,
                                      batch_size=batch_size, sampler=train_sampler, num_workers=num_workers,
                                      pin_memory=True)
        trainloader_iter = iter(trainloader)
        train_remain_sampler = data.sampler.SubsetRandomSampler(train_ids[partial_size:])
        trainloader_remain = data.DataLoader(train_dataset,
                                             batch_size=batch_size, sampler=train_remain_sampler, num_workers=1,
                                             pin_memory=True)
        trainloader_remain_iter = iter(trainloader_remain)
    else:
        train_sampler = data.sampler.SubsetRandomSampler(train_ids[:])
        trainloader = data.DataLoader(train_dataset,
                                      batch_size=batch_size, sampler=train_sampler, num_workers=num_workers,
                                      pin_memory=True)
        trainloader_iter = iter(trainloader)

    # Optimizer for segmentation network
    optimizer = build_optimizer(config['SOLVER']['OPTIMIZER'], model)
    optimizer.zero_grad()
    start_iteration = 0
    # resume & checkpoint
    if config['resume']:
        print("config['MODEL']['WEIGHTS'] is :",config['MODEL']['WEIGHTS'])
        checkpointer = DetectionCheckpointer(model, optimizer=optimizer)
        start_iteration = checkpointer.resume_or_load(config['MODEL']['WEIGHTS'], resume=True, train_again=config['train_again']).get("iteration", -1) + 1
    accumulated_loss_l = []
    if train_unlabeled:
        accumulated_loss_u = []
    if not os.path.exists(checkpoint_dir):
        os.makedirs(checkpoint_dir)
    with open(checkpoint_dir + '/config.json', 'w') as handle:
        json.dump(config, handle, indent=4, sort_keys=False)
    pickle.dump(train_ids, open(os.path.join(checkpoint_dir, 'train_split.pkl'), 'wb'))

    # epoch
    epochs_since_start = 0
    for i_iter in range(start_iteration, num_iterations):
        model.train()
        loss_l_value = 0
        if train_unlabeled:
            loss_u_value = 0
        optimizer.zero_grad()
        if lr_schedule:
            adjust_learning_rate(optimizer, i_iter)
        try:
            batch = next(trainloader_iter)
            if batch[0].shape[0] != batch_size:
                batch = next(trainloader_iter)
        except:
            epochs_since_start = epochs_since_start + 1
            print('Epochs since start: ', epochs_since_start)
            trainloader_iter = iter(trainloader)
            batch = next(trainloader_iter)
        weak_parameters = {"flip": 0}
        images, labels, _, _, _ = batch

        image_array = images.numpy()
        print(image_array.shape)  # (512, 512, 3)
        plt.imshow(image_array)  # 显示图片
        plt.axis('off')  # 不显示坐标轴
        plt.show()

        images = images.cuda()
        labels = labels.cuda()
        # images, labels = weakTransform(weak_parameters, data=images, target=labels)
        iter_rate = i_iter / num_iterations
        pred,gate_ori, gate_weight_ori = model(images, iter_rate)
        L_l = loss_calc(pred, labels)
        print('L_l:',L_l)
        # strong
        try:
            batch_remain = next(trainloader_remain_iter)
            if batch_remain[0].shape[0] != batch_size:
                batch_remain = next(trainloader_remain_iter)
        except:
            trainloader_remain_iter = iter(trainloader_remain)
            batch_remain = next(trainloader_remain_iter)
        images_remain, _, _, _, _ = batch_remain
        images_remain = images_remain.cuda()
        # images_remain = images
        inputs_u_w, _ = weakTransform(weak_parameters, data=images_remain)
        logits_u_w, gate_weak, gate_weight_weak = model(inputs_u_w, iter_rate)
        logits_u_w, _ = weakTransform(getWeakInverseTransformParameters(weak_parameters), data=logits_u_w.detach())
        softmax_u_w = torch.softmax(logits_u_w.detach(), dim=1)  # [1,19,768,768]
        max_probs, argmax_u_w = torch.max(softmax_u_w, dim=1)

        if mix_mask == "class":
            for image_i in range(batch_size):
                classes = torch.unique(argmax_u_w[image_i])
                classes = classes[classes != ignore_label]
                nclasses = classes.shape[0]
                classes = (classes[torch.Tensor(
                    np.random.choice(nclasses, int((nclasses - nclasses % 2) / 2), replace=False)).long()]).cuda()
                if image_i == 0:
                    MixMask = transformmasks.generate_class_mask(argmax_u_w[image_i], classes).unsqueeze(0).cuda()
                else:
                    MixMask = torch.cat((MixMask,
                                         transformmasks.generate_class_mask(argmax_u_w[image_i], classes).unsqueeze(
                                             0).cuda()))
        elif mix_mask == 'cut':
            img_size = inputs_u_w.shape[2:4]
            for image_i in range(batch_size):
                if image_i == 0:
                    MixMask = torch.from_numpy(transformmasks.generate_cutout_mask(img_size)).unsqueeze(
                        0).cuda().float()
                else:
                    MixMask = torch.cat((MixMask,
                                         torch.from_numpy(transformmasks.generate_cutout_mask(img_size)).unsqueeze(
                                             0).cuda().float()))
        elif mix_mask == "cow":
            img_size = inputs_u_w.shape[2:4]
            sigma_min = 8
            sigma_max = 32
            p_min = 0.5
            p_max = 0.5
            for image_i in range(batch_size):
                sigma = np.exp(np.random.uniform(np.log(sigma_min), np.log(sigma_max)))  # Random sigma
                p = np.random.uniform(p_min, p_max)  # Random p
                if image_i == 0:
                    MixMask = torch.from_numpy(
                        transformmasks.generate_cow_mask(img_size, sigma, p, seed=None)).unsqueeze(0).cuda().float()
                else:
                    MixMask = torch.cat((MixMask, torch.from_numpy(
                        transformmasks.generate_cow_mask(img_size, sigma, p, seed=None)).unsqueeze(
                        0).cuda().float()))
        elif mix_mask == None:
            MixMask = torch.ones((inputs_u_w.shape))
        strong_parameters = {"Mix": MixMask}
        if random_flip:
            strong_parameters["flip"] = random.randint(0, 1)
        else:
            strong_parameters["flip"] = 0
        if color_jitter:
            strong_parameters["ColorJitter"] = random.uniform(0, 1)
        else:
            strong_parameters["ColorJitter"] = 0
        if gaussian_blur:
            strong_parameters["GaussianBlur"] = random.uniform(0, 1)
        else:
            strong_parameters["GaussianBlur"] = 0

        inputs_u_s, _ = strongTransform(strong_parameters, data=images_remain)
        logits_u_s, gate_strong, gate_weight_strong = model(inputs_u_s, iter_rate)
        softmax_u_w_mixed, _ = strongTransform(strong_parameters, data=softmax_u_w)
        max_probs, pseudo_label = torch.max(softmax_u_w_mixed, dim=1)

        gate_weight_weak_diff = []
        gate_weight_strong_diff = []
        for i in range(len(gate_weight_ori)):
            for j in range(len(gate_weight_ori[i])):
                for k in range(len(gate_weight_ori[i][j])):
                    a = gate_weight_weak[i][j][k]
                    if gate_weight_ori[i][j][k] == []:
                        gate_weight_weak_diff.append(-1)
                        gate_weight_strong_diff.append(-1)
                    else:
                        gate_weight_weak_diff.append(round(gate_weight_weak[i][j][k][0] - gate_weight_ori[i][j][k][0],5))
                        gate_weight_strong_diff.append(round(gate_weight_strong[i][j][k][0] - gate_weight_ori[i][j][k][0],5))
        gate_weight_strong_diff_sum = np.sum(np.abs(np.array(gate_weight_strong_diff)))
        gate_weight_weak_diff_sum = np.sum(np.abs(np.array(gate_weight_weak_diff)))

        gate_weak_diff = []
        gate_strong_diff = []
        gate_weak_strong_diff = []
        for i in range(len(gate_ori)):
            for j in range(len(gate_ori[i])):
                    if gate_weak[i][j] is not gate_ori[i][j]:
                        gate_weak_diff.append(1)
                    else:
                        gate_weak_diff.append(0)
                    if gate_strong[i][j] is not gate_ori[i][j]:
                        gate_strong_diff.append(1)
                    else :
                        gate_strong_diff.append(0)
                    if gate_strong[i][j] is not gate_weak[i][j]:
                        gate_weak_strong_diff.append(1)
                    else :
                        gate_weak_strong_diff.append(0)
        gate_weak_diff_sum = np.sum(gate_weak_diff)
        gate_strong_diff_sum = np.sum(gate_strong_diff)
        gate_weak_strong_diff_sum = np.sum(gate_weak_strong_diff)

        print("finish")

    '''
        if pixel_weight == "threshold_uniform":
            unlabeled_weight = torch.sum(max_probs.ge(0.968).long() == 1).item() / np.size(
                np.array(pseudo_label.cpu()))
            pixelWiseWeight = unlabeled_weight * torch.ones(max_probs.shape).cuda()
        elif pixel_weight == "threshold":
            pixelWiseWeight = max_probs.ge(0.968).long().cuda()
        elif pixel_weight == 'sigmoid':
            max_iter = 10000
            pixelWiseWeight = sigmoid_ramp_up(i_iter, max_iter) * torch.ones(max_probs.shape).cuda()
        elif pixel_weight == False:
            pixelWiseWeight = torch.ones(max_probs.shape).cuda()
        # loss
        if consistency_loss == 'CE':
            L_u = consistency_weight * unlabeled_loss(logits_u_s, pseudo_label, pixelWiseWeight)
        elif consistency_loss == 'MSE':
            unlabeled_weight = torch.sum(max_probs.ge(0.968).long() == 1).item() / np.size(np.array(pseudo_label.cpu()))
            # softmax_u_w_mixed = torch.cat((softmax_u_w_mixed[1].unsqueeze(0),softmax_u_w_mixed[0].unsqueeze(0)))
            L_u = consistency_weight * unlabeled_weight * unlabeled_loss(logits_u_s, softmax_u_w_mixed)
        loss = L_l + L_u
    else:
        loss = L_l
    if len(gpus) > 1:
        loss = loss.mean()
        loss_l_value += L_l.mean().item()
        if train_unlabeled:
            loss_u_value += L_u.mean().item()
    else:
        loss_l_value += L_l.item()
        if train_unlabeled:
            loss_u_value += L_u.item()
    print("loss_u_value:",loss_u_value)
    '''

if __name__ == '__main__':

    print('---------------------------------Starting---------------------------------')
    args = get_arguments()
    config = json.load(open(args.config))
    print(args.config)
    gpus = (0, 1, 2, 3,4,5,6,7)[:args.gpus]
    model = config['model']
    dataset = config['dataset']

    if dataset == 'cityscapes':
        IMG_MEAN = np.flip(np.array([73.15835921, 82.90891754, 72.39239876]))
        num_classes = 19
        if config['training']['data']['split_id_list'] == 0:
            split_id = './splits/city/split_0.pkl'
        elif config['training']['data']['split_id_list'] == 1:
            split_id = './splits/city/split_1.pkl'
        elif config['training']['data']['split_id_list'] == 2:
            split_id = './splits/city/split_2.pkl'
        else:
            split_id = None

    elif dataset == 'pascal_voc':
        IMG_MEAN = np.flip(np.array([104.00698793, 116.66876762, 122.67891434]))
        num_classes = 21
        data_dir = './data/voc_dataset/'
        data_list_path = './data/voc_list/train_aug.txt'
        if config['training']['data']['split_id_list'] == 0:
            split_id = './splits/voc/split_0.pkl'
        else:
            split_id = None

    if config['pretrained'] == 'coco':
        restore_from = 'http://vllab1.ucmerced.edu/~whung/adv-semi-seg/resnet101COCO-41f33a49.pth'

    batch_size = config['training']['batch_size']
    num_iterations = config['training']['num_iterations']

    input_size_string = config['training']['data']['input_size']
    h, w = map(int, input_size_string.split(','))
    input_size = (h, w)

    ignore_label = config['ignore_label']  # 255 for PASCAL-VOC / 250 for Cityscapes

    learning_rate = config['training']['learning_rate']

    optimizer_type = config['training']['optimizer']
    lr_schedule = config['training']['lr_schedule']
    lr_power = config['training']['lr_schedule_power']
    weight_decay = config['training']['weight_decay']
    momentum = config['training']['momentum']
    num_workers = config['training']['num_workers']
    use_sync_batchnorm = config['training']['use_sync_batchnorm']
    random_seed = config['seed']

    labeled_samples = config['training']['data']['labeled_samples']

    # unlabeled CONFIGURATIONS
    train_unlabeled = config['train_unlabeled']
    mix_mask = config['training']['unlabeled']['mix_mask']
    pixel_weight = config['training']['unlabeled']['pixel_weight']
    consistency_loss = config['training']['unlabeled']['consistency_loss']
    consistency_weight = config['training']['unlabeled']['consistency_weight']
    random_flip = config['training']['unlabeled']['flip']
    color_jitter = config['training']['unlabeled']['color_jitter']
    gaussian_blur = config['training']['unlabeled']['blur']

    random_scale = config['training']['data']['scale']
    random_crop = config['training']['data']['crop']

    save_checkpoint_every = config['utils']['save_checkpoint_every']
    if args.resume:
        checkpoint_dir = os.path.join(*args.resume.split('/')[:-1]) + '_resume-' + start_writeable
    else:
        checkpoint_dir = os.path.join(config['utils']['checkpoint_dir'], start_writeable + '-' + args.name)

    log_dir = checkpoint_dir

    val_per_iter = config['utils']['val_per_iter']
    use_tensorboard = False # config['utils']['tensorboard']
    log_per_iter = config['utils']['log_per_iter']

    save_best_model = config['utils']['save_best_model']
    loss_weight = config['training']['unlabeled']['consistency_weight']

    if args.save_images:
        print('Saving unlabeled images')
        save_unlabeled_images = True
    else:
        save_unlabeled_images = False

    if config['model'] == 'dynamic_routing' and args.gpus == 1:
        config['MODEL']['BACKBONE']["NORM"] = 'BN'
        config['MODEL']['SEM_SEG_HEAD']["NORM"] = 'BN'

    main(config)


text:
from math import sin, cos, pi, sqrt, e, log, isclose, exp
from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors, cm, pyplot as plt
from scipy.constants import epsilon_0, mu_0
import scipy.constants
from fdtd_PCB_extensions.tissue import cole_cole_4, complex_permittivity_to_er_and_sigma, electric_field_penetration_depth
from mpl_toolkits.mplot3d import Axes3D
from propagation_lib import propagator

#old, deprecated, see _numerical_optimize.

#
# waveform = np.genfromtxt("NLTL_input_tone/NLTL-6275_cleaned.csv", delimiter=",", dtype=np.float, encoding='ascii', skip_header=0)
#
# waveform[:,0] = waveform[:,0] * 1e-12 #picoseonds
#
# input_frequency = 10e9
# sampling_frequency = 1000.0*input_frequency
#
# times = np.linspace(0, 1.0/input_frequency, 1000)
# input_data = np.interp(times, waveform[:,0], waveform[:,1])
#
# input_data -= np.mean(input_data) #remove the stupid zero reject
# input_data /= np.max(input_data) # normalize
#
# input_cycles = 10
#
# input_data = np.tile(input_data, 2)
# input_data[len(input_data/2):] = 0
#
# input_data = np.tile(input_data, input_cycles)
#

input_frequency = 9000.0e6
sampling_frequency = 1000.0*input_frequency
frequency_scale_factor = 1
times = np.linspace(0, 100*1.0/input_frequency, 100*1000)
# times = np.tile(times, 5)
input_data = np.zeros_like(times, dtype=np.float64)
input_data = np.sin(2.0*pi*times*input_frequency)
# input_data = np.tile(input_data, 2)

input_cycles = 2

p = propagator()
p.fourier_transform(input_data, sampling_frequency, 1.0)

depth_in_tissue = 0.1

p.populate_tissue_properties(48) #muscle

#
# for idx, depth_in_tissue in enumerate(np.arange(0, 0.1, 0.001)):
#     p.fourier_transform(input_data, sampling_frequency, 1.0)
#
#     spatial_phase = p.wavenumbers() * depth_in_tissue
#     # spatial_phase = np.zeros_like(p.mode_frequencies,dtype=np.float64)
#     p.attenuate(depth_in_tissue) # this global doesn't make any sense. pretty stupid design.
#
#     output = p.fourier_sum(spatial_phase)
#     output -= np.mean(output)
#     # print(depth_in_tissue, 1/np.max(output))
#     plt.plot(output[0:int(len(output)/(input_cycles/2))])
# plt.show()

# spatial_phase = p.wavenumbers() * depth_in_tissue
spatial_phase = np.zeros_like(p.mode_frequencies,dtype=np.float64)
p.attenuate(depth_in_tissue) # this global doesn't make any sense. pretty stupid design.

output = p.fourier_sum(spatial_phase)
output -= np.mean(output)
# print(depth_in_tissue, 1/np.max(output))
plt.plot(output)
plt.show()

#assert(output - input == 0)

# plt.plot(p.times, output)
#        plt.plot(p.times, input_data)

#Damping ratio: dimensionless number
#Damping constant:

osc_freq = 1000.0

# k = 1e-7
# charge = 1e7 * k

#turn the E field into a force
# p.a_n_coefficients *= charge
# p.b_n_coefficients *= charge

Q_factor = 100.0
# damping_force_coefficient = #N s/m

# mass = 120 * 1.66e-21 #kg
effective_mass = 120 * 1.66e-21 #kg #fix!

# beta = damping_force_coefficient/(2*effective_mass)
# beta = (osc_freq * 2*pi) / (2.0*Q_factor)
beta = 0.00000000000000001
# beta =/

# analytic_oscillator_scaling

print(p.analytic_oscillator_scaling(0.1, 1.0))

#beta is an exponential time constant, units /s.

print(p.mode_frequencies)
scale_factor = p.analytic_oscillator_scaling(beta, osc_freq)
print(scale_factor)

# plt.plot(p.mode_frequencies, p.a_n_coefficients)
plt.plot(p.mode_frequencies, scale_factor*p.a_n_coefficients)
plt.show()
#
p.a_n_coefficients *= scale_factor
p.b_n_coefficients *= scale_factor
osc_phase = p.analytic_oscillator_phase(beta, osc_freq)
oscillator_output = p.fourier_sum(spatial_phase + osc_freq)
#
# # oscillator_output /= 1e-12
#
plt.plot(p.times, oscillator_output)
plt.show()

reason: import numpy,import scipy,from scipy


text:
#!/usr/bin/env python
# Copyright 2017 Calico LLC
#
# 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.
# =========================================================================

from __future__ import print_function
from optparse import OptionParser
import os
import time

import h5py

import numpy as np
import pandas as pd
from scipy.stats import spearmanr
import tensorflow as tf

import matplotlib
matplotlib.use('PDF')
import matplotlib.pyplot as plt
import seaborn as sns

import basenji
"""basenji_predict_var.py

Make predictions and assess variance.
"""

################################################################################
# main
################################################################################
def main():
  usage = 'usage: %prog [options] <params_file> <model_file> <test_hdf5_file>'
  parser = OptionParser(usage)
  parser.add_option(
      '-d',
      dest='down_sample',
      default=1,
      type='int',
      help=
      'Down sample test computation by taking uniformly spaced positions [Default: %default]'
  )
  parser.add_option(
      '-m',
      dest='mc_n',
      default=0,
      type='int',
      help='Monte carlo iterations [Default: %default]')
  parser.add_option(
      '-o',
      dest='out_dir',
      default='predict_var',
      help='Output directory for test statistics [Default: %default]')
  parser.add_option(
      '-s',
      dest='save',
      default=False,
      action='store_true',
      help='Save predictions and variance arrays [Default: %default]')
  parser.add_option(
      '--rc',
      dest='rc',
      default=False,
      action='store_true',
      help=
      'Average the forward and reverse complement predictions when testing [Default: %default]'
  )
  parser.add_option(
      '-t',
      dest='targets',
      default=None,
      help='Comma-separated list of target indexes to plot [Default: %default]')
  parser.add_option(
      '-v',
      dest='valid',
      default=False,
      action='store_true',
      help='Process the validation set [Default: %default]')
  (options, args) = parser.parse_args()

  if len(args) != 3:
    parser.error('Must provide parameters, model, and test data HDF5')
  else:
    params_file = args[0]
    model_file = args[1]
    test_hdf5_file = args[2]

  if not os.path.isdir(options.out_dir):
    os.mkdir(options.out_dir)

  #######################################################
  # load data
  #######################################################
  data_open = h5py.File(test_hdf5_file)

  if not options.valid:
    test_seqs = data_open['test_in']
    test_targets = data_open['test_out']
    test_na = None
    if 'test_na' in data_open:
      test_na = data_open['test_na']

  else:
    test_seqs = data_open['valid_in']
    test_targets = data_open['valid_out']
    test_na = None
    if 'test_na' in data_open:
      test_na = data_open['valid_na']

  target_labels = [tl.decode('UTF-8') for tl in data_open['target_labels']]

  # limit targets
  if options.targets is None:
    target_indexes = np.arange(len(target_labels))
  else:
    target_indexes = np.array([int(ti) for ti in options.targets.split(',')])

  #######################################################
  # model parameters and placeholders

  job = basenji.dna_io.read_job_params(params_file)

  job['seq_length'] = test_seqs.shape[1]
  job['seq_depth'] = test_seqs.shape[2]
  job['num_targets'] = test_targets.shape[2]
  job['target_pool'] = int(np.array(data_open.get('pool_width', 1)))

  t0 = time.time()
  model = basenji.seqnn.SeqNN()
  model.build(job)
  print('Model building time %ds' % (time.time() - t0))

  # adjust for fourier
  job['fourier'] = 'train_out_imag' in data_open
  if job['fourier']:
    test_targets_imag = data_open['test_out_imag']
    if options.valid:
      test_targets_imag = data_open['valid_out_imag']

  #######################################################
  # predict

  # initialize batcher
  if job['fourier']:
    batcher_test = basenji.batcher.BatcherF(test_seqs, test_targets,
                                            test_targets_imag, test_na,
                                            model.batch_size, model.target_pool)
  else:
    batcher_test = basenji.batcher.Batcher(test_seqs, test_targets, test_na,
                                           model.batch_size, model.target_pool)

  # initialize saver
  saver = tf.train.Saver()

  with tf.Session() as sess:
    # load variables into session
    saver.restore(sess, model_file)

    # predict
    preds, preds_var = model.predict(
        sess,
        batcher_test,
        rc=options.rc,
        shifts=[-1, 0, 1],
        mc_n=options.mc_n,
        return_var=True,
        target_indexes=target_indexes,
        down_sample=options.down_sample)
    # preds, preds_var, preds_all = model.predict(sess, batcher_test, rc=options.rc, mc_n=options.mc_n, return_var=True, return_all=True, target_indexes=target_indexes, down_sample=options.down_sample)

    if options.save:
      np.save('%s/preds.npy' % options.out_dir, preds)
      np.save('%s/preds_var.npy' % options.out_dir, preds_var)
      # np.save('%s/preds_all.npy'%options.out_dir, preds_all)

  #######################################################
  # plot

  # plot prediction versus variance
  for tii in range(len(target_indexes)):
    preds_ti = preds[:, :, tii].flatten().astype('float32')
    preds_var_ti = preds_var[:, :, tii].flatten().astype('float32')
    out_pdf = '%s/t%d.pdf' % (options.out_dir, tii)

    cor, p = spearmanr(preds_ti, preds_var_ti)
    print('%4d  %.4f  %.4f' % (target_indexes[tii], cor, p))

    basenji.plots.jointplot(
        preds_ti,
        preds_var_ti,
        out_pdf,
        x_label='Prediction',
        y_label='Variance',
        figsize=(7, 7),
        sample=2000)

################################################################################
# __main__
################################################################################
if __name__ == '__main__':
  main()

max_stars_repo_name: egilbertson-ucsf/basenji

max_issues_repo_name: vaalessi/basenji

max_forks_repo_name: vaalessi/basenji


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

"""Tests for `magnipy` package."""

import pytest
import numpy as np
from click.testing import CliRunner
import matplotlib.pyplot as plt

@pytest.fixture
def response():
    """Sample pytest fixture.

    See more at: http://doc.pytest.org/en/latest/fixture.html
    """
    # import requests
    # return requests.get('https://github.com/audreyr/cookiecutter-pypackage')

def test_content(response):
    """Sample pytest test function with the pytest fixture as an argument."""
    # from bs4 import BeautifulSoup
    # assert 'GitHub' in BeautifulSoup(response.content).title.string

#def test_command_line_interface():
#    """Test the CLI."""
 #   runner = CliRunner()
#    result = runner.invoke(cli.main)
#    assert result.exit_code == 0
#    assert 'magnipy.cli.main' in result.output
#    help_result = runner.invoke(cli.main, ['--help'])
#    assert help_result.exit_code == 0
#    assert '--help  Show this message and exit.' in help_result.output

def test_profiles_nfw(plot=True):

    #x, y = np.linspace(-0.5, .5, 200), np.linspace(-0.5, 0.5, 200)

    x = np.loadtxt('xvalues.txt')
    y = np.linspace(0,0,len(x))
    #xx, yy = np.meshgrid(x, y)

    from MagniPy.MassModels.NFW import NFW
    from lenstronomy.LensModel.Profiles.nfw import NFW as NFW_L
    from MagniPy.MassModels.nfwT_temp import NFWt
    from MagniPy.MassModels.TNFW import TNFW
    import matplotlib.pyplot as plt

    rt = 1000
    fname = 'nfwdef_rt1000.txt'

    nfw = NFW()
    nfwL = NFW_L()
    TNFW = TNFW()
    nfwTL = NFWt()

    # nfw_params,nfw_params_L = nfw.params(x=0,y=0,mass=mass,mhm=mhm)
    # nfwt_params,nfwt_params_L = nfwT.params(x=0,y=0,mass=mass,mhm=mhm,trunc=rt)

    ks, rs = 0.1, 0.1
    theta_rs = 4 * ks * rs * (1 + np.log(.5))

    xdef, ydef = nfw.def_angle(x, y, center_x=0, center_y=0, theta_Rs=theta_rs, Rs=rs)
    xdef_t, ydef_t = TNFW.def_angle(x, y, center_x=0, center_y=0, theta_Rs=theta_rs, Rs=rs, r_trunc=rt)
    xdef_L, ydef_L = nfwL.derivatives(x=x, y=y, Rs=rs, theta_Rs=theta_rs)
    xdef_Lt, ydef_Lt = nfwTL.derivatives(x=x, y=y, Rs=rs, theta_Rs=theta_rs, t=rt)

    gravlens_xdef = np.loadtxt(fname)
    # xdef,ydef = nfw.def_angle(xx,yy,0,0,rs,ks)
    # xdef_L,ydef_L = nfwL.derivatives(x=xx,y=yy,**nfw_params_L)

    if plot:
        colors=['r','k','b','g']
        labels = ['my nfw','lenstronomy nfw','my nfw rt=1000rs','lensmodel']
        angles = [xdef,xdef_L,xdef_t,gravlens_xdef]
        for i,deflection in enumerate(angles[:-1]):
            plt.plot(x,deflection,color=colors[i],label=labels[i],alpha=0.5)
            plt.scatter(x, deflection, color=colors[i], label=labels[i], alpha=0.5)
        plt.plot(x,gravlens_xdef,color=colors[-1],label=labels[-1],alpha=0.5)
        plt.scatter(x, gravlens_xdef, color=colors[-1], label=labels[-1], alpha=0.5)
        plt.legend()
        plt.show()

    else:
        np.testing.assert_almost_equal(xdef, xdef_t, decimal=8)
        np.testing.assert_almost_equal(xdef, xdef_Lt, decimal=8)
        np.testing.assert_almost_equal(xdef, gravlens_xdef, decimal=8)
        np.testing.assert_almost_equal(xdef, xdef_L, decimal=8)

def test_profiles_nfwT(plot=True):

    rt = 0.4

    x = np.loadtxt('xvalues.txt')
    y = np.linspace(0, 0, len(x))

    from MagniPy.MassModels.NFW import NFW
    from lenstronomy.LensModel.Profiles.nfw import NFW as NFW_L
    from MagniPy.MassModels.nfwT_temp import NFWt
    from MagniPy.MassModels.TNFW import TNFW
    import matplotlib.pyplot as plt

    fname = 'nfwdef_rt.4.txt'

    nfw = NFW()
    nfwL = NFW_L()
    nfwT = TNFW()
    nfwTL = NFWt()

    # nfw_params,nfw_params_L = nfw.params(x=0,y=0,mass=mass,mhm=mhm)
    # nfwt_params,nfwt_params_L = nfwT.params(x=0,y=0,mass=mass,mhm=mhm,trunc=rt)

    ks, rs = 0.1, 0.1
    theta_rs = 4 * ks * rs * (1 + np.log(.5))

    xdef, ydef = nfw.def_angle(x, y, theta_Rs=theta_rs,Rs=rs,center_x=0, center_y=0)
    xdef_t, ydef_t = nfwT.def_angle(x, y, theta_Rs=theta_rs, Rs=rs, center_x=0, center_y=0, r_trunc=rt)
    xdef_L, ydef_L = nfwL.derivatives(x=x, y=y, Rs=rs, theta_Rs=theta_rs)
    xdef_Lt, ydef_Lt = nfwTL.derivatives(x=x, y=y, Rs=rs, theta_Rs=theta_rs, t=rt)
    gravlens_xdef = np.loadtxt(fname)
    # xdef,ydef = nfw.def_angle(xx,yy,0,0,rs,ks)
    # xdef_L,ydef_L = nfwL.derivatives(x=xx,y=yy,**nfw_params_L)

    if plot:
        colors=['r','k','b']
        labels = ['my nfw rt = 4*rs','lenstronomy nfw rt=4*rs','lensmodel']
        angles = [xdef_t,xdef_Lt,gravlens_xdef]
        for i,deflection in enumerate(angles[:-1]):
            plt.plot(x,deflection,color=colors[i],label=labels[i],alpha=0.5)
            plt.scatter(x, deflection, color=colors[i], label=labels[i], alpha=0.5)
        plt.plot(x,gravlens_xdef,color=colors[-1],label=labels[-1],alpha=0.5)
        plt.scatter(x, gravlens_xdef, color=colors[-1], label=labels[-1], alpha=0.5)
        plt.legend()
        plt.show()

    else:
        np.testing.assert_almost_equal(xdef_t, xdef_Lt, decimal=8)
        np.testing.assert_almost_equal(xdef_t, gravlens_xdef, decimal=8)

def test_solving_leq_single_plane():

    multiplane = False
    identifier = 'test_magnipy'

    from MagniPy.magnipy import Magnipy
    from MagniPy.lensdata import Data
    from MagniPy.LensBuild.main_deflector import Deflector
    from MagniPy.MassModels import SIE
    from MagniPy.Solver.solveroutines import SolveRoutines
    import matplotlib.pyplot as plt

    zmain,zsrc = 0.5,1.5
    srcx,srcy = -0.059,0.0065
    pos_sigma = [[0.003] * 4, [0.003] * 4]
    flux_sigma = [.4] * 4
    tdel_sigma = np.array([0, 2, 2, 2]) * 10
    sigmas = [pos_sigma, flux_sigma, tdel_sigma]
    solver = SolveRoutines(0.5, 1,5)
    init = Magnipy(0.5, 1.5, use_lenstronomy_halos=True)

    datax = [-0.6894179,0.4614137,0.05046926,0.3233231]
    datay = [-0.1448493,0.3844455,0.5901034,-0.4361114]
    datam = [ 0.42598153,1,0.75500025,0.36081329]
    datat = [ 0,6.314742,6.777318,9.38729]

    start = {'R_ein': .7, 'ellip': .05, 'ellip_theta': 0, 'x': 0, 'y': 0, 'shear': 0, 'shear_theta': 0}
    start = Deflector(subclass=SIE.SIE(), redshift=zmain, tovary=True,
                      varyflags=['1', '1', '1', '1', '1', '0', '0', '0', '0', '0'], **start)

    data_to_fit = [Data(x=datax,y=datay,m=datam,t=datat,source=[srcx,srcy])]

    profiles = ['NFW']

    for profile in profiles:

        halo_models2 = ["plaw_"+profile+"_[.005]_6_10_[0,1]_['uniformnfw',[3,500]]_1"]

        realizationsdata = init.generate_halos(halo_models2,Nrealizations=1)

        dx,dy = -.01,0
        if profile=='NFW':
            ks = 0.01
            newparams = {'x':data_to_fit[0].x[1]+dx,'y':data_to_fit[0].y[1]+dy,'ks':ks,'rs':0.02}

        realizationsdata[0][0].update(method='lensmodel',**newparams)
        realizationsdata[0] = [realizationsdata[0][0]]

        optimized_lenstronomy,systems = solver.fit_src_plane(macromodel=start, datatofit=data_to_fit[0],
                                                             realizations=realizationsdata, multiplane=multiplane,
                                                             method='lenstronomy', ray_trace=True, sigmas=sigmas,
                                                             identifier=identifier, srcx=srcx, srcy=srcy, grid_rmax=.05,
                                                             res=.001, source_shape='GAUSSIAN', source_size=0.0012)

        optimized_lensmodel,systems = solver.two_step_optimize(macromodel=start, datatofit=data_to_fit[0],
                                                               realizations=realizationsdata, multiplane=multiplane,
                                                               method='lensmodel', ray_trace=True, sigmas=sigmas,
                                                               identifier=identifier, srcx=srcx, srcy=srcy,
                                                               grid_rmax=.05,
                                                               res=.001, source_shape='GAUSSIAN', source_size=0.0012)

        np.testing.assert_almost_equal(optimized_lensmodel[0].m, optimized_lenstronomy[0].m, decimal=2)
        np.testing.assert_almost_equal(optimized_lensmodel[0].x, optimized_lenstronomy[0].x, decimal=4)
        np.testing.assert_almost_equal(optimized_lensmodel[0].y, optimized_lenstronomy[0].y, decimal=4)
        np.testing.assert_almost_equal(optimized_lensmodel[0].x, optimized_lenstronomy[0].x, decimal=4)
        #plt.scatter(optimized_lenstronomy[0].x,optimized_lenstronomy[0].y,color='r',marker='x',alpha=0.5)
        #lt.scatter(optimized_lensmodel[0].x,optimized_lensmodel[0].y,color='k',alpha=0.5)
        #plt.scatter(data_to_fit[0].x,data_to_fit[0].y,color='k',marker='+',s=70,alpha=0.5)
        #plt.scatter(newparams['x'],newparams['y'],color='m',s=50)
        print optimized_lensmodel[0].m
        print optimized_lenstronomy[0].m
        #plt.show()

def test_SIE():
    from MagniPy.MassModels.SIE import SIE

    values = np.loadtxt('../../SIE_def.txt',skiprows=5)
    sie = SIE()

    xdef = values[:,3]
    ydef = values[:,4]
    x = values[:,0]
    y = values[:,1]

    theta_E = 0.7
    ellip = 0.45
    ellip_theta = -10

    center_x,center_y = 0,0
    plt.plot(x,xdef,color='k')

    siexdef = sie.def_angle(x,y,center_x=center_x,center_y=center_y,theta_E=theta_E,q=1-ellip,phi_G=ellip_theta*np.pi*180**-1)[0]
    plt.plot(x,siexdef,color='r')
    plt.show()

def test_Shear():
    from MagniPy.MassModels.ExternalShear import Shear
    s = Shear()

    shear,shear_theta = 0.02,40

    values = np.loadtxt('../../shear_def.txt',skiprows=5)
    xdef = values[:, 3]
    ydef = values[:, 4]
    x = values[:, 0]
    y = values[:, 1]
    plt.plot(x,xdef)
    plt.plot(x,s.def_angle(x,y,shear,shear_theta)[0])

    plt.show()

def test_SIEShear():
    from MagniPy.MassModels.SIE import SIE
    from MagniPy.MassModels.ExternalShear import Shear
    s = Shear()

    shear, shear_theta = 0.02, 40

    values_shear = np.loadtxt('../../shear_def.txt', skiprows=5)

    x = values_shear[:, 0]
    y = values_shear[:, 1]

    values_sie = np.loadtxt('../../SIE_def.txt', skiprows=5)
    sie = SIE()

    values_sieshear = np.loadtxt('../../sieshear_def.txt',skiprows=5)

    x = values_sie[:, 0]
    y = values_sie[:, 1]

    theta_E = 0.7
    ellip = 0.45
    ellip_theta = -10

    center_x, center_y = 0, 0
    plt.plot(x, values_sieshear[:,3], color='k')

    siexdef = sie.def_angle(x, y, center_x=center_x, center_y=center_y, theta_E=theta_E, q=1 - ellip,
                            phi_G=ellip_theta * np.pi * 180 ** -1)[0]
    shearxdef = s.def_angle(x,y,shear,shear_theta)[0]
    plt.plot(x,siexdef+shearxdef,color='r')
    plt.show()

def test_nfws():
    from lenstronomy.LensModel.Profiles.nfw import NFW
    x_loc = np.linspace(0,.1,100)
    y_loc = np.linspace(0,0,100)
    args = {'theta_Rs':0.01,'Rs':0.02,'center_x':0,'center_y':0.03}
    nfw = NFW()
    xdefL, ydefL = nfw.derivatives(x=x_loc, y=y_loc, **args)
    from MagniPy.MassModels.NFW import NFW
    nfw = NFW()
    xplus,yplus = nfw.def_angle(x=x_loc,y=y_loc,**args)
    print ydefL
    print yplus
    np.testing.assert_almost_equal(ydefL, yplus, decimal=5)

def los_test():
    from MagniPy.Solver.GravlensWrap.generate_input import Header
    from MagniPy.Solver.GravlensWrap._call import call_lensmodel
    from MagniPy.util import polar_to_cart

    h = Header(zlens=0.5,zsrc=1.5,omega_M=0.3,hval=0.7)
    lines = h.inputstring
    e1,e2 = polar_to_cart(0.2,-80)
    sievals = 'alpha 1 0 0 '+str(e1)+' '+str(e2)+' '+'0 0 0 0 1 0.5\n'
    nfwvals = ' nfw '+'0.050387342258 0.92 0.4 0 0 0 0 1.39523097511 0 0 0.2\n'
    lines += 'setlens 2 1 1\n '+sievals+nfwvals+'0 0 0 0 0 0 0 0 0 0\n0 0 0 0 0 0 0 0 0 0\n'
    outfile='losnfw_def.txt'
    lines += 'plotdef1 '+outfile+' 0.8 1.2 1000 0.4 0.4 1000'

    with open('../../los_nfw.txt','w') as f:
        f.write(lines)
    #call_lensmodel('los_nfw.txt','../../')

    values = np.loadtxt('../../losnfw_def.txt',skiprows=5)
    x,xdef_lensmodel,ydef_lensmodel = values[:,0],values[:,3],values[:,4]
    plt.plot(x,xdef_lensmodel,color='k')
    plt.show()

#los_test()
#test_nfws()
#test_SIEShear()
#test_profiles_nfwT(plot=True)
#test_profiles_nfw(plot=True)
#test_solving_leq_single_plane()

max_stars_repo_name: dangilman/MagniPy

max_issues_repo_name: dangilman/MagniPy

max_forks_repo_name: dangilman/MagniPy


text:
import numpy as np
import quatMath
import compFilterFuncs
from scipy.spatial.transform import Rotation
from scipy.spatial.transform import Slerp
from scipy.spatial.transform import RotationSpline
from scipy import signal
from numpy.linalg import inv

class complementaryFilter:

    def __init__(self, kGravity, gamma_, kAngularVelocityThreshold, kAccelerationThreshold,
               kDeltaAngularVelocityThreshold, gain_acc_, bias_alpha_, do_bias_estimation_,
               do_adaptive_gain_,initialized_,steady_state_,timeStep,time, sensorType):
    # define CONSTANT params
        self.kGravity = kGravity
        self.gamma_ = gamma_
        self.kAngularVelocityThreshold = kAngularVelocityThreshold
        self.kAccelerationThreshold = kAccelerationThreshold
        self.kDeltaAngularVelocityThreshold = kDeltaAngularVelocityThreshold
        self.gain_acc_ =gain_acc_
        self.bias_alpha_ = bias_alpha_
        self.do_bias_estimation_ = do_bias_estimation_
        self.do_adaptive_gain_ = do_adaptive_gain_
        self.initialized_ = initialized_
        self.steady_state_ = steady_state_
        self.timeStep = timeStep
        self.time = time
        self.sensorType = sensorType #if zero: Xsens, if 1: HL

        self.q0_ = 1     #these can be changed with xsens quaternions
        self.q1_ = 0
        self.q2_ = 0
        self.q3_ = 0

        self.wx_prev_ = 0
        self.wy_prev_ = 0
        self.wz_prev_ = 0
        self.wx_bias_ = 0
        self.wy_bias_ = 0
        self.wz_bias_ = 0

        self.new_wx = np.zeros(len(time), dtype=float)
        self.new_wy = np.zeros(len(time), dtype=float)
        self.new_wz = np.zeros(len(time), dtype=float)
        self.neworienta= np.zeros(len(time), dtype=float)
        self.neworientb= np.zeros(len(time), dtype=float)
        self.neworientc= np.zeros(len(time), dtype=float)
        self.neworientd= np.zeros(len(time), dtype=float)

        self.gx  = np.zeros(len(time), dtype = float)
        self.gy  = np.zeros(len(time), dtype = float)
        self.gz  = np.zeros(len(time), dtype = float)

    def execute(self,accX, accY, accZ, gyrX, gyrY,gyrZ,
                interp_hlAccX, interp_hlAccY, interp_hlAccZ, interp_hlGyrX, interp_hlGyrY, interp_hlGyrZ):

        if self.sensorType == 0:
            iteratorType = accX
        elif self.sensorType ==1:
            iteratorType = interp_hlAccX

        for i in range(len(iteratorType)):
            #x,y,z are world coords, z points DOWN, x pointing EAST, y pointing NORTH (ENU frame)
            # if self.sensorType == 0:  #xsens
            #     ax = accX[i]
            #     ay = accY[i]
            #     az = accZ[i]

            #     wx = gyrX[i]
            #     wy = gyrY[i]
            #     wz = gyrZ[i]

            # if self.sensorType ==1 : #hl
            #     ax = interp_hlAccX[i]
            #     ay = interp_hlAccY[i]
            #     az = -interp_hlAccZ[i]

            #     wx = interp_hlGyrX[i]
            #     wy = interp_hlGyrY[i]
            #     wz = interp_hlGyrZ[i]

            if self.sensorType == 0:  #xsens
                ax = accZ[i]
                ay = -accY[i]
                az = accX[i]

                wx = gyrZ[i]
                wy = -gyrY[i]
                wz = gyrX[i]

            if self.sensorType ==1 : #hl
                ax = interp_hlAccZ[i]
                ay = -interp_hlAccY[i]
                az = interp_hlAccX[i]

                wx = interp_hlGyrZ[i]
                wy = -interp_hlGyrY[i]
                wz = interp_hlGyrX[i]

            if (self.initialized_ == False):
                # start up
                self.q0_,self.q1_,self.q2_,self.q3_ = compFilterFuncs.getMeasurement(ax, ay, az,                                                                 self.q0_, self.q1_, self.q2_, self.q3_)
                self.initialized_ = True

            self.wx_bias_, self.wy_bias_, self.wz_bias_, self.wx_prev_, self.wy_prev_, self.wz_prev_ = compFilterFuncs.updateBiases(ax,ay,az,wx,wy,wz,
                                                                    self.kAccelerationThreshold,self.kAngularVelocityThreshold,
                                                                    self.kDeltaAngularVelocityThreshold,self.kGravity,
                                                                    self.wx_bias_, self.wy_bias_, self.wz_bias_,
                                                                    self.wx_prev_, self.wy_prev_, self.wz_prev_,
                                                                    self.bias_alpha_)

            q0_pred = 0
            q1_pred = 0
            q2_pred = 0
            q3_pred = 0
            q0_pred, q1_pred, q2_pred, q3_pred = compFilterFuncs.Prediction(wx, wy, wz, self.timeStep, q0_pred, q1_pred, q2_pred, q3_pred,
                                                                     self.q0_, self.q1_, self.q2_, self.q3_,
                                                                     self.wx_bias_, self.wy_bias_, self.wz_bias_)

            dq0_acc = 0
            dq1_acc = 0
            dq2_acc = 0
            dq3_acc = 0

            dq0_acc,dq1_acc,dq2_acc,dq3_acc,self.gx[i], self.gy[i], self.gz[i] = compFilterFuncs.AccelComp(ax, ay, az,
                                                     q0_pred, q1_pred, q2_pred, q3_pred,
                                                     dq0_acc, dq1_acc, dq2_acc, dq3_acc)

            #adaptive gain constant
            gain = self.gain_acc_

            #FINALIZE RESULTS
            dq0_acc, dq1_acc, dq2_acc, dq3_acc=  quatMath.scaleQuaternion(gain, dq0_acc, dq1_acc, dq2_acc, dq3_acc)

            self.q0_, self.q1_, self.q2_, self.q3_ =   quatMath.quaternionMultiplication(q0_pred, q1_pred, q2_pred, q3_pred,
                                                                    dq0_acc, dq1_acc, dq2_acc, dq3_acc,
                                                                    self.q0_, self.q1_, self.q2_, self.q3_);

            self.q0_, self.q1_, self.q2_, self.q3_ =   quatMath.normalizeQuaternion(self.q0_, self.q1_, self.q2_, self.q3_)

            self.neworienta[i] = self.q0_
            self.neworientb[i] = self.q1_
            self.neworientc[i] = self.q2_
            self.neworientd[i]= self.q3_

            ##handle angular velocity bias

            self.new_wx[i] = wx - self.wx_bias_
            self.new_wy[i] = wy - self.wy_bias_
            self.new_wz[i] = wz - self.wz_bias_

        return (self.neworienta, self.neworientb, self.neworientc, self.neworientd,
                self.new_wx, self.new_wy, self.new_wz,self.gx, self.gy, self.gz)

    def solveRotationMatrix(self,neworienta, neworientb, neworientc, neworientd,time):
        stackedMatrix = np.zeros((len(time)+1,3,3))
        for i in range(len(time)):
            r = Rotation.from_quat([neworienta[i],neworientb[i],neworientc[i],neworientd[i]])
            rotmat = r.as_matrix()
            stackedMatrix[i,:,:]= rotmat

        return(stackedMatrix)

    def solveOrientation(self,xsensstackedMatrix,hlstackedMatrix,xsensFPS):

        #jump the first five seconds because the person needs to
        #start the streams and get into a stable position
        #after the first 5 seconds, an averaging of the first second
        #with euler angles can be obtained and from this a rotation matrix can be created
        jump = int(np.floor(60* xsensFPS))
        XS_eulX=np.zeros(int(xsensFPS), dtype=np.float64)
        XS_eulY=np.zeros(int(xsensFPS), dtype=np.float64)
        XS_eulZ=np.zeros(int(xsensFPS), dtype=np.float64)
        HL_eulX=np.zeros(int(xsensFPS), dtype=np.float64)
        HL_eulY=np.zeros(int(xsensFPS), dtype=np.float64)
        HL_eulZ=np.zeros(int(xsensFPS), dtype=np.float64)
        eulerX=np.zeros(int(xsensFPS), dtype=np.float64)
        eulerY=np.zeros(int(xsensFPS), dtype=np.float64)
        eulerZ=np.zeros(int(xsensFPS), dtype=np.float64)

        #from xsens to world, and then from world to hl
        R_hl_xsens = np.zeros((len(xsensstackedMatrix)-1,3,3))
        for i in range(len(xsensstackedMatrix)-1):
            R_hl_xsens[i,:,:] = np.matmul(inv(hlstackedMatrix[i,:,:]),xsensstackedMatrix[i,:,:])

        for i in range(int(xsensFPS)):

            temp = Rotation.from_matrix(R_hl_xsens[i+jump])
            r = temp.as_euler('zyx', degrees=True)
            eulerX[i] = r[2]
            eulerY[i] = r[1]
            eulerZ[i] = r[0]

        for i in range(int(xsensFPS)):
            temp = Rotation.from_matrix(xsensstackedMatrix[i+jump])
            r = temp.as_euler('zyx', degrees=True)
            XS_eulX[i] = r[2]
            XS_eulY[i] = r[1]
            XS_eulZ[i] = r[0]

        for i in range(int(xsensFPS)):
            temp = Rotation.from_matrix(hlstackedMatrix[i+jump])
            r = temp.as_euler('zyx', degrees=True)
            HL_eulX[i] = r[2]
            HL_eulY[i] = r[1]
            HL_eulZ[i] = r[0]

        #average the results to optimize the orientation between them
        avg_hl_eulX = np.mean(HL_eulX)
        avg_hl_eulY = np.mean(HL_eulY)
        avg_hl_eulZ = np.mean(HL_eulZ)
        avg_xs_eulX = np.mean(XS_eulX)
        avg_xs_eulY = np.mean(XS_eulY)
        avg_xs_eulZ = np.mean(XS_eulZ)

        avg_eulX = np.mean(eulerX)
        avg_eulY = np.mean(eulerY)
        avg_eulZ = np.mean(eulerZ)

        #create a new rotation matrix in between xsens and the hl

        #rot matrix that takes values from HL to world
        avg_R_hl = Rotation.from_euler('zyx',[avg_hl_eulZ, avg_hl_eulY, avg_hl_eulX], degrees=True)
        avg_R_hl = avg_R_hl.as_matrix()

        #rot matrix that takes values from XSENS to world
        avg_R_xs = Rotation.from_euler('zyx',[avg_xs_eulZ, avg_xs_eulY, avg_xs_eulX], degrees=True)
        avg_R_xs = avg_R_xs.as_matrix()

        #rot matrix that takes values from Xsens to HL

#        avg_R_hl_xsens = np.matmul(inv(avg_R_hl),avg_R_xs)

        avg_R_hl_xsens = Rotation.from_euler('zyx',[avg_eulX, avg_eulY, avg_eulZ], degrees = True)
        avg_R_hl_xsens = avg_R_hl_xsens.as_matrix()

        #kalibr
 #        avg_R_hl_xsens = np.array([
 #            [-0.95052692,  0.29522118, -0.09665936],
 # [-0.2675078,  -0.61972534,  0.73782118],
 # [ 0.15791818,  0.72717602,  -0.66803958]
 # ])

        return(R_hl_xsens, avg_R_hl_xsens)

    def XsensOnHL(self,coarseAligned_xsensGyrX, coarseAligned_xsensGyrY, coarseAligned_xsensGyrZ,
                  coarseAligned_xsensAccX, coarseAligned_xsensAccY, coarseAligned_xsensAccZ, avg_R_hl_xsens):
        '''
        transfers the xsens gyro values on HL for cross-correlation
        '''

        xsensStacked_omega = np.array([coarseAligned_xsensGyrX,
                                       coarseAligned_xsensGyrY,
                                       coarseAligned_xsensGyrZ
            ])
        xsensStacked_alpha = np.array([coarseAligned_xsensAccX,
                                       coarseAligned_xsensAccY,
                                       coarseAligned_xsensAccZ
            ])
        #copy the original stack to a var
        xsensStacked_omega_copy = xsensStacked_omega
        xsensStacked_alpha_copy = xsensStacked_alpha

        xsens_Stacked_omega = np.matmul(avg_R_hl_xsens, xsensStacked_omega)
        coarseAligned_xsensGyrX = xsens_Stacked_omega[0]
        coarseAligned_xsensGyrY = xsens_Stacked_omega[1]
        coarseAligned_xsensGyrZ = xsens_Stacked_omega[2]

        xsens_Stacked_alpha = np.matmul(avg_R_hl_xsens, xsensStacked_alpha)
        coarseAligned_xsensAccX = xsens_Stacked_alpha[0]
        coarseAligned_xsensAccY = xsens_Stacked_alpha[1]
        coarseAligned_xsensAccZ = xsens_Stacked_alpha[2]

        return (coarseAligned_xsensGyrX, coarseAligned_xsensGyrY, coarseAligned_xsensGyrZ,
                coarseAligned_xsensAccX, coarseAligned_xsensAccY, coarseAligned_xsensAccZ,
                xsensStacked_omega_copy,xsensStacked_alpha_copy)

max_stars_repo_name: ozgunkaratas/mastersthesis

max_issues_repo_name: ozgunkaratas/mastersthesis

max_forks_repo_name: ozgunkaratas/mastersthesis

reason: import numpy,from numpy,from scipy


text:
#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
@author: Cem Rıfkı Aydın
"""
import numpy as np
import spacy
from GloVe.glove_generator import generate_glove_vecs

from SVD_U.svd_U_generator import generate_svd_u_vecs

nlp = spacy.load("en_core_web_sm")

def generate_average_word_vecs_from_corpus(corp, word_vecs):
    """
    This function generates document embeddings in the corpus by averaging all the
    word vectors therein.

    :param corp: A list of corpus sentences that are composed of words.
    :type corp: list
    :param word_vecs: A dictionary containing words as keys and the corresponding vectors as values.
    :type word_vecs: dict
    :return: Averaged document vectors. This corresponds to all the corpus reviews.
    :rtype: np.array
    """

    mean_ = np.array(list(word_vecs.values())).mean(0)

    upd_corp = []

    for index, review in enumerate(corp):
        if index + 1 % 1000 == 0:
            print(index)
        mean_review_vec = generate_average_word_vec(review, word_vecs, mean_)
        upd_corp.append(mean_review_vec)
    return np.array(upd_corp)

def generate_average_word_vec(review, word_vecs, mean_):
    """
    This function generates a document embedding by averaging all the vectors of words
        occurring in that document.

    :param review: A text (e.g. review or e-mail) of str type.
    :type review: str
    :param word_vecs: A dictionary containing words as keys and the corresponding vectors as values.
    :type word_vecs: dict
    :param mean_: For the non-existent words, the mean vector of the whole training vectors is used.
    :type mean_: np.array or list
    :return: The averaged word vector.
    :rtype: np.array
    """

    review_vectors = np.array([word_vecs[tok] if tok in word_vecs else mean_ for tok in review])
    # review_vectors = np.array([word_vecs[tok.text] if tok.text in word_vecs else mean_ for tok in nlp(review)])
    mean_review_vec = review_vectors.mean(0)
    return mean_review_vec

def _similarity_query(word, number, word_vectors):
    """
    This function returns the most similar words to a given query.

    :param word: The target word.
    :type word: str
    :param number: The number of the most similar words to the above-given word.
    :type number: int
    :param word_vectors: word_vectors (): The set of word vectors, whereby the similarities between
            words are to be computed.
    :type word_vectors: dict
    :return: The most similar tokens to the given word.
    :rtype: list
    """

    keys, vals = list(word_vectors.keys()), list(word_vectors.values())

    word_vec = word_vectors[word]
    dst = (np.dot(vals, word_vec)
           / np.linalg.norm(vals, axis=1)
           / np.linalg.norm(word_vec))
    word_ids = np.argsort(-dst)

    return [keys[x] for x in word_ids[:number]
            if word in keys]

def generate_word_embeddings(embedding_type, revs):
    """
    This main function generates either GloVe or SVD - U embeddings.

    :param embedding_type: The type of vectors. This can be either glove or SVD_U.
    :type embedding_type: str
    :param revs: The list of reviews.
    :type revs: list or np.array
    :return: Word vectors modelled according to one of the algorithms stated above.
    :rtype: dict
    """

    import time
    start_time = time.time()

    type_dict = {"glove": generate_glove_vecs,
                 "SVD_U": generate_svd_u_vecs}

    word_vectors = type_dict[embedding_type](revs)

    elapsed_time = time.time() - start_time

    print("Word vectors have been generated. Elapsed time for it is shown below:")

    print(time.strftime("%H:%M:%S", time.gmtime(elapsed_time)))

    # print(_similarity_query("more", 5, word_vectors))
    # print(_similarity_query("cat", 5, word_vectors))
    # print(_similarity_query("best", 5, word_vectors))
    # print(_similarity_query("movie", 5, word_vectors))

    return word_vectors


text:
import os
import sys
import tensorflow as tf
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib import rcParams
import numpy as np
import seaborn as sns
from keras.models import load_model
from sklearn.preprocessing import StandardScaler

sns.set(style='whitegrid', palette='muted', font_scale=1.5)

rcParams['figure.figsize'] = 22, 10

if len(sys.argv) < 7:
    print("Too few arguments")
    quit()
elif len(sys.argv) > 7:
    print("Too many arguments")
    quit()

# Read the arguments
for i in range (1, len(sys.argv)):
    if(sys.argv[i] == "-d"):
        dataset_name = sys.argv[i+1]
    elif(sys.argv[i] == "-n"):
        num_time_series = int(sys.argv[i+1])
    elif(sys.argv[i] == "-mae"):
        mae = float(sys.argv[i+1])

num = 100
csv_path = os.path.join(os.path.abspath(__file__), "../../../dir/")
csv_path = os.path.join(csv_path, dataset_name)

model_path = os.path.join(os.path.abspath(__file__), "../../../models/detect_model.h5")

df = pd.read_csv(csv_path, header=None, delimiter='\t')
file_ids = df.iloc[:, [0]].values
df = df.drop(df.columns[0], axis=1)
df = df.transpose()

TIME_STEPS = 30

# Training data: 95%, Test data: 5%
train_size = int(len(df) * 0.95)
test_size = len(df) - train_size

scaler = StandardScaler()

# Creating a data structure with 30 time-steps and 1 output
X_train = []
y_train = []
for step in range (0, num):

    training_set = df.iloc[:train_size, [step]].values
    training_set_scaled = scaler.fit_transform(training_set)

    for i in range(TIME_STEPS, train_size):
        X_train.append(training_set_scaled[i-TIME_STEPS:i, 0])
        y_train.append(training_set_scaled[i, 0])

X_train, y_train = np.array(X_train), np.array(y_train)
X_train = np.reshape(X_train, (X_train.shape[0], X_train.shape[1], 1))

dataset_test = df.iloc[train_size:, [0]]

# Load model
model = load_model(model_path)

X_train_pred = model.predict(X_train)
train_mae_loss = np.mean(np.abs(X_train_pred - X_train), axis=1)

for series in range(0, num_time_series):
    dataset_train = df.iloc[:train_size, [series]]
    dataset_test = df.iloc[train_size:, [series]]
    dataset_total = pd.concat((dataset_train, dataset_test), axis = 0)

    inputs = dataset_total[len(dataset_total) - len(dataset_test) - TIME_STEPS:].values
    inputs = inputs.reshape(-1,1)
    inputs = scaler.transform(inputs)
    X_test = []
    for i in range(TIME_STEPS, test_size):
        X_test.append(inputs[i-TIME_STEPS:i, 0])

    X_test = np.array(X_test)
    X_test = np.reshape(X_test, (X_test.shape[0], X_test.shape[1], 1))

    X_test_pred = model.predict(X_test)
    test_mae_loss = np.mean(np.abs(X_test_pred - X_test), axis=1)

    THRESHOLD = mae

    # DataFrame containing the loss and the anomalies (values above the threshold)
    test_score_df = pd.DataFrame(index=dataset_test[TIME_STEPS:].index)
    test_score_df['loss'] = test_mae_loss
    test_score_df['threshold'] = THRESHOLD
    test_score_df['anomaly'] = test_score_df.loss > test_score_df.threshold
    test_score_df['value'] = dataset_test[TIME_STEPS:]

    anomalies = test_score_df[test_score_df.anomaly == True]

    plt.plot(
        dataset_test[TIME_STEPS:].index,
        dataset_test[TIME_STEPS:],
        label='value')

    sns.scatterplot(
        anomalies.index,
        anomalies.value,
        color=sns.color_palette()[3],
        s=52,
        label='anomaly')

    plt.xticks(rotation=25)
    plt.legend()
    fig = plt.gcf()
    fig.set_size_inches(18.5, 10.5)
    plt.savefig(os.path.join(os.path.abspath(__file__), "../../../dir/exports/q2/graph"+ str(series) +".png"))
    plt.clf()