text:
from io import BytesIO
import time
import picamera
import numpy as np
import BlackWhite2
from PIL import Image

output = 'output.jpg'

def getArray():
    with picamera.PiCamera() as camera:
       camera.framerate = 24
       time.sleep(2)
       camera.capture(output)
       print('picture taken')
       outputnew = BlackWhite2.binarize_image(output, 75)
       print(outputnew)
       return outputnew

if __name__ == "__main__":
    getArray()


text:
# coding: UTF-8
# adapted from https://github.com/arosh/BM25Transformer/blob/master/bm25.py

import numpy as np
import scipy as sp

from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.utils.validation import check_is_fitted
from sklearn.feature_extraction.text import _document_frequency
from utils import *

class BM25Transformer(BaseEstimator, TransformerMixin):
    """
    Parameters
    ----------
    use_idf : boolean, optional (default=True)
    k1 : float, optional (default=2.0)
    b  : float, optional (default=0.75)
    References
    ----------
    Okapi BM25: a non-binary model - Introduction to Information Retrieval
    http://nlp.stanford.edu/IR-book/html/htmledition/okapi-bm25-a-non-binary-model-1.html
    """
    def __init__(self, use_idf=True, k1=2.0, b=0.75):
        self.use_idf = use_idf
        self.k1 = k1
        self.b = b

    def fit(self, X):
        """
        Parameters
        ----------
        X : sparse matrix, [n_samples, n_features] document-term matrix
        """
        if not sp.sparse.issparse(X):
            X = sp.sparse.csc_matrix(X)
        if self.use_idf:
            n_samples, n_features = X.shape
            df = _document_frequency(X)
            idf = np.log((n_samples - df + 0.5) / (df + 0.5))
            self._idf_diag = sp.sparse.spdiags(idf, diags=0, m=n_features, n=n_features)

        doc_len = X.sum(axis=1)
        self._average_document_len = np.average(doc_len)

        return self

    def transform(self, X, copy=True):
        """
        Parameters
        ----------
        X : sparse matrix, [n_samples, n_features] document-term matrix
        copy : boolean, optional (default=True)
        """
        if hasattr(X, 'dtype') and np.issubdtype(X.dtype, np.float):
            # preserve float family dtype
            X = sp.sparse.csr_matrix(X, copy=copy)
        else:
            # convert counts or binary occurrences to floats
            X = sp.sparse.csr_matrix(X, dtype=np.float, copy=copy)

        n_samples, n_features = X.shape

        # Document length (number of terms) in each row
        # Shape is (n_samples, 1)
        doc_len = X.sum(axis=1)
        # Number of non-zero elements in each row
        # Shape is (n_samples, )
        sz = X.indptr[1:] - X.indptr[0:-1]

        # In each row, repeat `doc_len` for `sz` times
        # Shape is (sum(sz), )
        # Example
        # -------
        # dl = [4, 5, 6]
        # sz = [1, 2, 3]
        # rep = [4, 5, 5, 6, 6, 6]
        rep = np.repeat(np.asarray(doc_len), sz)

        # Compute BM25 score only for non-zero elements
        nom = self.k1 + 1
        denom = X.data + self.k1 * (1 - self.b + self.b * rep / self._average_document_len)
        data = X.data * nom / denom

        X = sp.sparse.csr_matrix((data, X.indices, X.indptr), shape=X.shape)

        if self.use_idf:
            check_is_fitted(self, '_idf_diag', 'idf vector is not fitted')

            expected_n_features = self._idf_diag.shape[0]
            if n_features != expected_n_features:
                raise ValueError("Input has n_features=%d while the model"
                                 " has been trained with n_features=%d" % (
                                     n_features, expected_n_features))
            X = X * self._idf_diag

        return X

with timer('merge additional files'):
    train = merge_breed_name(train)

orig_cols = train.columns
train[text_features].fillna('missing', inplace=True)
with timer('tfidf_bm25 + svd / nmf'):
    vectorizer = make_pipeline(
        TfidfVectorizer(),
        BM25Transformer(use_idf=True, k1=2.0, b=0.75),
        TruncatedSVD(n_components=n_components, random_state=seed)
    )
    X = vectorizer.fit_transform(train['Description'])
    X = pd.DataFrame(X, columns=['tfidf_bm25_svd_{}'.format(i) for i in range(n_components)])
    train = pd.concat([train, X], axis=1)
    del vectorizer;
    gc.collect()

with timer('count_bm25 + svd / nmf'):
    vectorizer = make_pipeline(
        CountVectorizer(),
        BM25Transformer(use_idf=True, k1=2.0, b=0.75),
        TruncatedSVD(n_components=n_components, random_state=seed),
    )
    X = vectorizer.fit_transform(train['Description'])
    X = pd.DataFrame(X, columns=['count_bm25_svd_{}'.format(i) for i in range(n_components)])
    train = pd.concat([train, X], axis=1)
    del vectorizer;
    gc.collect()

with timer('tfidf2_bm25 + svd / nmf'):
    vectorizer = make_pipeline(
        TfidfVectorizer(min_df=2, max_features=20000,
                        strip_accents='unicode', analyzer='word', token_pattern=r'\w{1,}',
                        ngram_range=(1, 3), use_idf=1, smooth_idf=1, sublinear_tf=1, stop_words='english'),
        BM25Transformer(use_idf=True, k1=2.0, b=0.75),
        TruncatedSVD(n_components=n_components, random_state=seed)
    )
    X = vectorizer.fit_transform(train['Description'])
    X = pd.DataFrame(X, columns=['tfidf2_bm25_svd_{}'.format(i) for i in range(n_components)])
    train = pd.concat([train, X], axis=1)
    del vectorizer;
    gc.collect()

with timer('count2_bm25 + svd / nmf'):
    vectorizer = make_pipeline(
        CountVectorizer(min_df=2, max_features=20000,
                        strip_accents='unicode', analyzer='word', token_pattern=r'\w{1,}',
                        ngram_range=(1, 3), stop_words='english'),
        BM25Transformer(use_idf=True, k1=2.0, b=0.75),
        TruncatedSVD(n_components=n_components, random_state=seed)
    )
    X = vectorizer.fit_transform(train['Description'])
    X = pd.DataFrame(X, columns=['count2_bm25_svd_{}'.format(i) for i in range(n_components)])
    train = pd.concat([train, X], axis=1)
    del vectorizer;
    gc.collect()
new_cols = [c for c in train.columns if c not in orig_cols]
train[new_cols].to_feather("../feature/bm25.feather")

ext: py

reason: import numpy,import scipy


text:
import numba as _numba
import numpy as _np
import numexpr as _ne
from ..util import rotation as _rotation

@_numba.njit(parallel=True)
def _decaysum(a, t, tau):
    # reference without scan
    x = _np.zeros_like(a)
    for j in _numba.prange(x.shape[0]):
        x[j, 0] = a[j, 0]
        for i in range(1, x.shape[1]):
            decay = _np.exp(-(t[j, i] - t[j, i - 1]) / tau)
            x[j, i] = a[j, i] + decay * x[j, i - 1]
    return x

@_numba.njit(parallel=True)
def _psort(a, b):
    '''
    idx = np.argsort(a, axis=-1)
    a = a[np.arange(a.shape[0])[:, None], idx]
    b = b[np.arange(a.shape[0])[:, None], idx]
    '''
    for i in _numba.prange(len(a)):
        ids = _np.argsort(a[i])
        a[i, :] = a[i][ids]
        b[i, :] = b[i][ids]

def _ab2(x):
    return x.imag ** 2 + x.real ** 2

def _integral(amp, t0, tau):
    """
    integrates the exponential decay of |x|^2 with supports t0 and decay tau
    """
    t0s = _np.atleast_2d(t0)
    amps = _np.atleast_2d(amp)
    _psort(t0s, amps)
    i = _ab2(_decaysum(amps, t0s, tau))
    td = _np.diff(t0s, axis=-1)
    return _np.sum(-tau / 2 * i[:, :-1] * _np.expm1(-2 * td / tau), axis=-1) + tau / 2 * i[:, -1]

def simulate(simobject, Ndet, pixelsize, detz, k, c, tau, pulsewidth, detangles=(0, 0), settings='', threads=None, *args, **kwargs):
    """
    Time dependent simulation with decaying amplitudes (cpu version with few optimisations).
    simobject: simobject to use for simulation (in lengthunit)
    pixelsize: pixelsize (in lengthunit)
    detz: Detector-sample distance
    k: angular wave number (in 1/lengthunit)
    c: speed of light in (lengthunit/timeunit)
    tau: decay time (in timeunit)
    pulsewidth: FWHM of gaussian exciation pulse (in timeunit)
    settings: string, can contain
        scale - do 1/r intensity scaling
    threads: blocksize to do multiple pixels at once (memory tradeoff)
    """

    from math import sin, cos

    if _np.size(Ndet) == 1:
        Ndet = [Ndet, Ndet]
    if 'scale' in settings:
        eq = 'exp(-1j*(s*k-phases))/d'
    else:
        eq = 'exp(-1j*(s*k-phases))'
    n = _np.prod(Ndet)
    if threads is None:
        threads = 4
    # do highest power of two of n <= threads pixels at once. tradeoff between memory allocations and call overhead.
    blocksize = min(threads, (n & (~(n - 1))))

    dets = _np.array(
        _np.meshgrid(pixelsize * (_np.arange(Ndet[0]) - (Ndet[0] / 2)), pixelsize * (_np.arange(Ndet[1]) - (Ndet[1] / 2)), detz)
    ).T.reshape(-1, 3)
    if _np.any(detangles):
        dets = dets @ _rotation(*detangles, 0).T

    res = _np.zeros(Ndet[0] * Ndet[1]).reshape(-1, blocksize)
    data = simobject.get()
    times = _np.random.randn(simobject.N) * (pulsewidth / 2.35)

    # this could be parallelised but array sharing might be annoying to implement, so either mpi run everything or be happy with the threaded parts
    for j, det in enumerate(dets.reshape(-1, blocksize, 3)):
        d = _np.linalg.norm(det, axis=-1)[:, None]
        q = det / d
        q[..., -1] -= 1  # incoming beam along z direction
        phases = data[:, 3]
        s = _np.inner(q, data[:, :3])  # path difference
        e = _ne.evaluate(eq)  # complex e field
        t = -s / c + times.T  # arrival time
        res[j] = _integral(e, t, tau)
    res = res.reshape(Ndet[0], Ndet[1])
    return res

def simulate_gen(
    simobject, Ndet, pixelsize, detz, k, c, tau, pulsewidth, detangles=(0, 0), settings='', maximg=_np.inf, threads=None, *args, **kwargs
):
    """
    Time dependent simulation with decaying amplitudes (cpu version with few optimisations).
    simobject: simobject to use for simulation (in lengthunit)
    pixelsize: pixelsize (in lengthunit)
    detz: Detector-sample distance
    k: angular wave number (in 1/lengthunit)
    c: speed of light in (lengthunit/timeunit)
    tau: decay time (in timeunit)
    pulsewidth: FWHM of gaussian exciation pulse (in timeunit)
    settings: string, can contain
        scale - do 1/r intensity scaling
    threads: blocksize to do multiple pixels at once (memory tradeoff)
    """

    # as each simulation is slow compared to setup, we dont currently do any work ahead and just call simulate()
    i = 0
    while i < maximg:
        yield simulate(simobject, Ndet, pixelsize, detz, k, c, tau, pulsewidth, detangles, settings, threads, *args, **kwargs)
        i += 1

reason: import numpy,import numba,import numexpr


text:
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 13 18:41:38 2021

@author: divyoj
"""
## importing libraries:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation
import os
# # note that this must be executed before 'import numba'
# os.environ['NUMBA_DISABLE_INTEL_SVML'] = '1'
from numba import njit
import time as process_time
import plotting
from scipy.integrate import solve_ivp
@njit
def do_timestep(t,z,aT,bT,alpha, beta, gamma, zeta):
    ''' function to give dxdt at a time step '''

    aL = z[0*(nx*ny):1*(nx*ny)].reshape((ny,nx))
    bL = z[1*(nx*ny):2*(nx*ny)].reshape((ny,nx))
    aR = z[2*(nx*ny):3*(nx*ny)].reshape((ny,nx))
    bR = z[3*(nx*ny):4*(nx*ny)].reshape((ny,nx))

    # total membrane protein concentration:
    a0 = aL + aR
    b0 = bL + bR
    # intitating dxdt to arrays of zeros:
    daL=np.zeros((ny,nx));daR=np.zeros((ny,nx));dbL=np.zeros((ny,nx));dbR=np.zeros((ny,nx));

    ## Equations for al,aR,bl and bR:
    daL=alpha*(aT-a0)-beta*aL+beta*gamma*( aL*np.roll(bR,+1) ) -zeta*(aL-aR)**3;
    # aR
    daR=alpha*(aT-a0)-beta*aR+beta*gamma*( aR*np.roll(bL,-1) ) -zeta*(aR-aL)**3;
    # bL
    dbL=alpha*(bT-b0)-beta*bL+beta*gamma*( bL*np.roll(aR,+1) ) -zeta*(bL-bR)**3;
    # bR
    dbR=alpha*(bT-b0)-beta*bR+beta*gamma*( bR*np.roll(aL,-1) ) -zeta*(bR-bL)**3;

    # # checking if membrane concentration went above total concetration:

    daL=daL*((aT>=a0) | (daL<0));  daR=daR*((aT>=a0) | (daR<0))
    dbL=dbL*((bT>=b0) | (dbL<0));  dbR=dbR*((bT>=b0) | (dbR<0))

    daL=daL*((aL>=0) | (daL>0));       daR=daR*((aR>=0) | (daR>0))
    dbL=dbL*((bL>=0) | (dbL>0));       dbR=dbR*((bR>=0) | (dbR>0))

    return np.concatenate((daL.flatten(),dbL.flatten(),daR.flatten(),dbR.flatten()))

#@njit
def simulate(rho,alpha, beta, gamma, zeta):
    ''' function to iterate over time and return arrays with the result '''
    ## initilizing the arrays to store the values over time:
    aL_t = np.zeros((T_max+1,ny,nx)); aR_t = np.zeros((T_max+1,ny,nx));
    bL_t = np.zeros((T_max+1,ny,nx)); bR_t = np.zeros((T_max+1,ny,nx));

    # total proteins in the cells
    aT = rho + np.zeros((ny,nx))
    bT = rho + np.zeros((ny,nx))

    ## initializing aL,bR,bL,aR
    aL = np.zeros((ny,nx)) + 0.1*rho; aR = np.zeros((ny,nx)) + 0.100001*rho
    bL = np.zeros((ny,nx)) + 0.100001*rho; bR = np.zeros((ny,nx)) + 0.1*rho

    ## Collecting the initial conditions into a single array:
    ic = np.array(aL.flatten().tolist()+bL.flatten().tolist()+aR.flatten().tolist()+bR.flatten().tolist())

    ## Solving the initial value problem:
    sol = solve_ivp(lambda t,y: do_timestep(t,y,aT,bT,alpha, beta, gamma, zeta),t_span=[0,T_max],y0=ic,method='RK45',t_eval=list(np.linspace(0,T_max,T_max+1)))

    t = sol.t
    aball = sol.y

    for t_index, ts in enumerate(t):

        aball_at_ts = aball[:,t_index]

        aL_t[t_index]= aball_at_ts[0*(nx*ny):1*(nx*ny)].reshape((ny,nx));
        bL_t[t_index]= aball_at_ts[1*(nx*ny):2*(nx*ny)].reshape((ny,nx));
        aR_t[t_index]= aball_at_ts[2*(nx*ny):3*(nx*ny)].reshape((ny,nx));
        bR_t[t_index]= aball_at_ts[3*(nx*ny):4*(nx*ny)].reshape((ny,nx));

    return (aL_t,bL_t,aR_t,bR_t)

if __name__ == "__main__":
    # Lattice:
    w,h = 50,2;
    dx,dy=0.01,0.1;
    nx=int(w/dx)
    ny=1;#int(h/dx);

    # time:
    T_max=200;

    # parameters:
    alpha=10;
    gamma=1;beta=1;
    zeta=0.1;

    #%% Characterisation over rho:
    alpha_array=[5,10,20]# np.round(np.logspace(-1,2,5),2)
    f, axs = plt.subplots(3,1,figsize=(4,9))

    for alphai, alpha in enumerate(alpha_array):
        print("alpha=",alpha)
        # folder for storing the data:
        folder="./Uniform_rho/"+"zeta="+str(zeta)+"_alpha="+str(alpha)+"/"
        if not os.path.exists(folder):
            os.makedirs(folder)

        rhos = np.around(np.linspace(0.8,1.5,50),3);rho0_array=rhos*1;
        p_array=rhos*1;delta_p_array=rhos*1

        for ri, rho in enumerate(rhos):
            print (ri, rho)

            aL_t, bL_t, aR_t, bR_t = simulate(rho,alpha, beta, gamma, zeta)

             # Plotting at each rho
            #%% Plotting at each rho:
            rho0_array[ri],p_array[ri],delta_p_array[ri]=plotting.plots_at_rho(aL_t,bL_t,aR_t,bR_t,rho,folder)

        ## rho0 vs rho
        axs[0].plot(rhos,rho0_array,label=str(alpha));
        axs[0].set_title(r"$\rho_{0} \ v/s \  \rho$");
        axs[0].set_ylabel(r"$\rho_{0}$")
        axs[0].set_xlabel(r"$\rho$")
        axs[0].legend(ncol=2)
        ## rho0 vs rho
        axs[1].plot(rhos,p_array,label=str(alpha));
        axs[1].set_title(r'$<p> \ v/s \ \rho$');
        axs[1].set_ylabel(r"$<p>$")
        axs[1].set_xlabel(r"$\rho$")
        axs[1].legend(ncol=2)
        ## rho0 vs rho
        axs[2].plot(rhos,delta_p_array,label=str(alpha));
        axs[2].set_title(r'$\Delta p \ v/s \ \rho $');
        axs[2].set_ylabel(r"$\Delta p $")
        axs[2].set_xlabel(r"$\rho$")
        axs[2].legend(ncol=2)

    f.suptitle(r"zeta="+str(zeta))
    f.subplots_adjust(top=0.85, bottom=0.20, left=0.20, right=0.95, hspace=0.50,wspace=0.50)

    f.savefig("./Uniform_rho/"+"Uniform_over_rho_zeta="+str(zeta)+".png",dpi=500)
    plt.close()

    #%% Characterisation over rho:
    rhos_array=[0.8,1,1.5]# np.round(np.logspace(-1,2,5),2)
    f, axs = plt.subplots(3,1,figsize=(4,9))

    for rhoi, rho in enumerate(rhos_array):
        print("rho=",rho)
        #folder for storing the data:
        folder="./Uniform_alpha/"+"zeta="+str(zeta)+"_rho="+str(rho)+"/"
        if not os.path.exists(folder):
            os.makedirs(folder)

        alphas = np.round(np.logspace(-1, 1.5, num=30),3);rho0_array=alphas.copy();
        p_array=alphas.copy();delta_p_array=alphas.copy()

        for alphai, alpha in enumerate(alphas):
            print (alphai, alpha)

            aL_t, aR_t, bL_t, bR_t = simulate(rho,alpha, beta, gamma, zeta)

            # Plotting at each rho
            rho0_array[alphai],p_array[alphai],delta_p_array[alphai]=plotting.plots_at_rho(aL_t,bL_t,aR_t,bR_t,alpha,folder)

            ## rho0 vs alpha
        axs[0].plot(alphas,rho0_array,label=str(rho));
        axs[0].set_title(r"$\rho_{0} \ vs\  \alpha$");
        axs[0].set_ylabel(r"$\rho_{0}$")
        axs[0].set_xscale('log')
        axs[0].legend(ncol=2)

        axs[1].plot(alphas,p_array,label=str(rho));
        axs[1].set_title(r'$<p> \ vs \  \alpha$');
        axs[1].set_ylabel(r"$<p>$")
        axs[1].set_xscale('log')
        axs[1].legend(ncol=2)

        axs[2].plot(alphas,delta_p_array,label=str(rho));
        axs[2].set_title(r'$\Delta p \ vs\  \alpha$');
        axs[2].set_ylabel(r"$\Delta p$")
        axs[2].set_xlabel(r"$\alpha$")
        axs[2].set_xscale('log')
        axs[2].legend(ncol=2)

    f.suptitle(r"zeta="+str(zeta))
    f.subplots_adjust(top=0.85, bottom=0.20, left=0.20, right=0.95, hspace=0.50,wspace=0.50)

    f.savefig("./Uniform_alpha/"+"Uniform_over_alpha_zeta="+str(zeta)+".png",dpi=500)
    plt.close()

reason: import numpy,from scipy,import numba,from numba


text:
'''
Test that maximizing likelihood allows to correctly recover distribution parameters for all
distributions exposed to the user.
'''
# Standard library imports
from typing import Iterable, List, Tuple

# Third-party imports
import mxnet as mx
import numpy as np
import pytest
from pydantic import PositiveFloat, PositiveInt

# First-party imports
from gluonts.model.common import NPArrayLike
from gluonts.distribution.box_cox_tranform import (
    InverseBoxCoxTransform,
    InverseBoxCoxTransformOutput,
)
from gluonts.distribution import (
    DistributionOutput,
    StudentT,
    StudentTOutput,
    MultivariateGaussian,
    MultivariateGaussianOutput,
    LowrankMultivariateGaussian,
    LowrankMultivariateGaussianOutput,
    NegativeBinomial,
    NegativeBinomialOutput,
    Laplace,
    LaplaceOutput,
    Gaussian,
    GaussianOutput,
    PiecewiseLinear,
    PiecewiseLinearOutput,
    Binned,
    BinnedOutput,
)
from gluonts.distribution.transformed_distribution_output import (
    TransformedDistributionOutput
)
from gluonts.distribution.transformed_distribution import (
    TransformedDistribution
)

NUM_SAMPLES = 2000
BATCH_SIZE = 32
TOL = 0.2
START_TOL_MULTIPLE = 1

np.random.seed(1)
mx.random.seed(1)

def inv_softplus(y: NPArrayLike) -> np.ndarray:
    # y = log(1 + exp(x))  ==>  x = log(exp(y) - 1)
    return np.log(np.exp(y) - 1)

def maximum_likelihood_estimate_sgd(
    distr_output: DistributionOutput,
    samples: mx.ndarray,
    init_biases: List[mx.ndarray.NDArray] = None,
    num_epochs: PositiveInt = 5,
    learning_rate: PositiveFloat = 1e-2,
    hybridize: bool = True,
) -> Iterable[float]:
    model_ctx = mx.cpu()

    arg_proj = distr_output.get_args_proj()
    arg_proj.initialize()

    if hybridize:
        arg_proj.hybridize()

    if init_biases is not None:
        for param, bias in zip(arg_proj.proj, init_biases):
            param.params[param.prefix + "bias"].initialize(
                mx.initializer.Constant(bias), force_reinit=True
            )

    trainer = mx.gluon.Trainer(
        arg_proj.collect_params(),
        'sgd',
        {'learning_rate': learning_rate, 'clip_gradient': 10.0},
    )

    dummy_data = mx.nd.array(np.ones((len(samples), 1)))

    train_data = mx.gluon.data.DataLoader(
        mx.gluon.data.ArrayDataset(dummy_data, samples),
        batch_size=BATCH_SIZE,
        shuffle=True,
    )

    dummy_out = mx.nd.array(np.ones((1, 1)))

    for e in range(num_epochs):
        cumulative_loss = 0
        num_batches = 0
        # inner loop
        for i, (data, sample_label) in enumerate(train_data):
            # data: (batch, 1), the "1" defines the dimension of the projection layer and should be kept
            data = data.as_in_context(model_ctx)
            # sample_label: (batch,), for univariate and (batch, d) for multivariate
            if sample_label.shape[-1] == 1:
                sample_label = sample_label.as_in_context(model_ctx).squeeze(
                    axis=-1
                )
            else:
                sample_label = sample_label.as_in_context(model_ctx)

            with mx.autograd.record():
                distr_args = arg_proj(data)
                distr = distr_output.distribution(distr_args)
                loss = distr.loss(sample_label)
            loss.backward()
            trainer.step(BATCH_SIZE)
            num_batches += 1

            cumulative_loss += mx.nd.mean(loss).asscalar()
        print("Epoch %s, loss: %s" % (e, cumulative_loss / num_batches))
    return [param[0].asnumpy() for param in arg_proj(dummy_out)]

@pytest.mark.parametrize("mu_sigma_nu", [(2.3, 0.7, 6.0)])
@pytest.mark.parametrize("hybridize", [True, False])
def test_studentT_likelihood(
    mu_sigma_nu: Tuple[float, float, float], hybridize: bool
) -> None:
    '''
    Test to check that maximizing the likelihood recovers the parameters
    '''
    # test instance
    mu, sigma, nu = mu_sigma_nu

    # generate samples
    mus = mx.nd.zeros((NUM_SAMPLES, 1)) + mu
    sigmas = mx.nd.zeros((NUM_SAMPLES, 1)) + sigma
    nus = mx.nd.zeros((NUM_SAMPLES, 1)) + nu

    distr = StudentT(mus, sigmas, nus)
    samples = distr.sample()

    # nu takes very long to learn, so we initialize it at the true value.
    # transform used is softplus(x) + 2
    init_bias = [
        mu - START_TOL_MULTIPLE * TOL * mu,
        inv_softplus(sigma - START_TOL_MULTIPLE * TOL * sigma),
        inv_softplus(nu - 2),
    ]

    mu_hat, sigma_hat, nu_hat = maximum_likelihood_estimate_sgd(
        StudentTOutput(),
        samples,
        init_biases=init_bias,
        hybridize=hybridize,
        num_epochs=PositiveInt(10),
        learning_rate=1e-2,
    )

    assert (
        np.abs(mu_hat - mu) < TOL * mu
    ), f"mu did not match: mu = {mu}, mu_hat = {mu_hat}"
    assert (
        np.abs(sigma_hat - sigma) < TOL * sigma
    ), f"sigma did not match: sigma = {sigma}, sigma_hat = {sigma_hat}"
    assert (
        np.abs(nu_hat - nu) < TOL * nu
    ), "nu0 did not match: nu0 = %s, nu_hat = %s" % (nu, nu_hat)

@pytest.mark.parametrize("mu_sigma", [(1.0, 0.1)])
@pytest.mark.parametrize("hybridize", [True, False])
def test_gaussian_likelihood(mu_sigma: Tuple[float, float], hybridize: bool):
    '''
    Test to check that maximizing the likelihood recovers the parameters
    '''
    # test instance
    mu, sigma = mu_sigma

    # generate samples
    mus = mx.nd.zeros((NUM_SAMPLES, 1)) + mu
    sigmas = mx.nd.zeros((NUM_SAMPLES, 1)) + sigma

    distr = Gaussian(mus, sigmas)
    samples = distr.sample()

    init_biases = [
        mu - START_TOL_MULTIPLE * TOL * mu,
        inv_softplus(sigma - START_TOL_MULTIPLE * TOL * sigma),
    ]

    mu_hat, sigma_hat = maximum_likelihood_estimate_sgd(
        GaussianOutput(),
        samples,
        init_biases=init_biases,
        hybridize=hybridize,
        learning_rate=PositiveFloat(0.001),
        num_epochs=PositiveInt(5),
    )

    assert (
        np.abs(mu_hat - mu) < TOL * mu
    ), f"mu did not match: mu = {mu}, mu_hat = {mu_hat}"
    assert (
        np.abs(sigma_hat - sigma) < TOL * sigma
    ), f"alpha did not match: sigma = {sigma}, sigma_hat = {sigma_hat}"

@pytest.mark.timeout(10)
def test_multivariate_gaussian() -> None:
    num_samples = 2000
    dim = 2

    mu = np.arange(0, dim) / float(dim)

    L_diag = np.ones((dim,))
    L_low = 0.1 * np.ones((dim, dim)) * np.tri(dim, k=-1)
    L = np.diag(L_diag) + L_low
    Sigma = L.dot(L.transpose())

    distr = MultivariateGaussian(mu=mx.nd.array(mu), L=mx.nd.array(L))

    samples = distr.sample(num_samples)

    mu_hat, L_hat = maximum_likelihood_estimate_sgd(
        MultivariateGaussianOutput(dim=dim),
        samples,
        init_biases=None,  # todo we would need to rework biases a bit to use it in the multivariate case
        hybridize=False,
        learning_rate=PositiveFloat(0.01),
        num_epochs=PositiveInt(10),
    )

    distr = MultivariateGaussian(
        mu=mx.nd.array([mu_hat]), L=mx.nd.array([L_hat])
    )

    Sigma_hat = distr.variance[0].asnumpy()

    assert np.allclose(
        mu_hat, mu, atol=0.1, rtol=0.1
    ), f"mu did not match: mu = {mu}, mu_hat = {mu_hat}"
    assert np.allclose(
        Sigma_hat, Sigma, atol=0.1, rtol=0.1
    ), f"Sigma did not match: sigma = {Sigma}, sigma_hat = {Sigma_hat}"

@pytest.mark.timeout(10)
def test_lowrank_multivariate_gaussian() -> None:
    num_samples = 2000
    dim = 2
    rank = 1

    mu = np.arange(0, dim) / float(dim)
    D = np.eye(dim) * (np.arange(dim) / dim + 0.5)
    W = np.sqrt(np.ones((dim, rank)) * 0.2)
    Sigma = D + W.dot(W.transpose())

    distr = LowrankMultivariateGaussian(
        mu=mx.nd.array([mu]),
        D=mx.nd.array([np.diag(D)]),
        W=mx.nd.array([W]),
        dim=dim,
        rank=rank,
    )

    assert np.allclose(
        distr.variance[0].asnumpy(), Sigma, atol=0.1, rtol=0.1
    ), f"did not match: sigma = {Sigma}, sigma_hat = {distr.variance[0]}"

    samples = distr.sample(num_samples).squeeze().asnumpy()

    mu_hat, D_hat, W_hat = maximum_likelihood_estimate_sgd(
        LowrankMultivariateGaussianOutput(dim=dim, rank=rank),
        samples,
        learning_rate=PositiveFloat(0.01),
        num_epochs=PositiveInt(10),
        init_biases=None,  # todo we would need to rework biases a bit to use it in the multivariate case
        hybridize=False,
    )

    distr = LowrankMultivariateGaussian(
        dim=dim,
        rank=rank,
        mu=mx.nd.array([mu_hat]),
        D=mx.nd.array([D_hat]),
        W=mx.nd.array([W_hat]),
    )

    Sigma_hat = distr.variance.asnumpy()

    assert np.allclose(
        mu_hat, mu, atol=0.2, rtol=0.1
    ), f"mu did not match: mu = {mu}, mu_hat = {mu_hat}"

    assert np.allclose(
        Sigma_hat, Sigma, atol=0.1, rtol=0.1
    ), f"alpha did not match: sigma = {Sigma}, sigma_hat = {Sigma_hat}"

@pytest.mark.parametrize("mu", [6.0])
@pytest.mark.parametrize("hybridize", [True, False])
def test_deterministic_l2(mu: float, hybridize: bool) -> None:
    '''
    Test to check that maximizing the likelihood recovers the parameters.
    This tests uses the Gaussian distribution with fixed variance and sample mean.
    This essentially reduces to determistic L2.
    '''
    # generate samples
    mu = mu
    mus = mx.nd.zeros(NUM_SAMPLES) + mu

    deterministic_distr = Gaussian(mu=mus, sigma=0.1 * mx.nd.ones_like(mus))
    samples = deterministic_distr.sample()

    class GaussianFixedVarianceOutput(GaussianOutput):
        @classmethod
        def domain_map(cls, F, mu, sigma):
            sigma = 0.1 * F.ones_like(sigma)
            return mu.squeeze(axis=-1), sigma.squeeze(axis=-1)

    mu_hat, _ = maximum_likelihood_estimate_sgd(
        GaussianFixedVarianceOutput(),
        samples,
        init_biases=[3 * mu, 0.1],
        hybridize=hybridize,
        num_epochs=PositiveInt(1),
    )

    assert (
        np.abs(mu_hat - mu) < TOL * mu
    ), f"mu did not match: mu = {mu}, mu_hat = {mu_hat}"

@pytest.mark.parametrize("mu", [1.0])
@pytest.mark.parametrize("hybridize", [True, False])
def test_deterministic_l1(mu: float, hybridize: bool) -> None:
    '''
    Test to check that maximizing the likelihood recovers the parameters.
    This tests uses the Laplace distribution with fixed variance and sample mean.
    This essentially reduces to determistic L1.
    '''
    # generate samples
    mu = mu
    mus = mx.nd.zeros(NUM_SAMPLES) + mu

    class LaplaceFixedVarianceOutput(LaplaceOutput):
        @classmethod
        def domain_map(cls, F, mu, b):
            b = 0.1 * F.ones_like(b)
            return mu.squeeze(axis=-1), b.squeeze(axis=-1)

    deterministic_distr = Laplace(mu=mus, b=0.1 * mx.nd.ones_like(mus))
    samples = deterministic_distr.sample()

    mu_hat, _ = maximum_likelihood_estimate_sgd(
        LaplaceFixedVarianceOutput(),
        samples,
        init_biases=[3 * mu, 0.1],
        learning_rate=1e-3,
        hybridize=hybridize,
    )

    assert (
        np.abs(mu_hat - mu) < TOL * mu
    ), f"mu did not match: mu = {mu}, mu_hat = {mu_hat}"

@pytest.mark.parametrize("mu_alpha", [(2.5, 0.7)])
@pytest.mark.parametrize("hybridize", [True, False])
def test_neg_binomial(mu_alpha: Tuple[float, float], hybridize: bool) -> None:
    '''
    Test to check that maximizing the likelihood recovers the parameters
    '''
    # test instance
    mu, alpha = mu_alpha

    # generate samples
    mus = mx.nd.zeros((NUM_SAMPLES, 1)) + mu
    alphas = mx.nd.zeros((NUM_SAMPLES, 1)) + alpha

    neg_bin_distr = NegativeBinomial(mu=mus, alpha=alphas)
    samples = neg_bin_distr.sample()

    init_biases = [
        inv_softplus(mu - START_TOL_MULTIPLE * TOL * mu),
        inv_softplus(alpha + START_TOL_MULTIPLE * TOL * alpha),
    ]

    mu_hat, alpha_hat = maximum_likelihood_estimate_sgd(
        NegativeBinomialOutput(),
        samples,
        hybridize=hybridize,
        init_biases=init_biases,
        num_epochs=PositiveInt(15),
    )

    assert (
        np.abs(mu_hat - mu) < TOL * mu
    ), f"mu did not match: mu = {mu}, mu_hat = {mu_hat}"
    assert (
        np.abs(alpha_hat - alpha) < TOL * alpha
    ), f"alpha did not match: alpha = {alpha}, alpha_hat = {alpha_hat}"

@pytest.mark.timeout(10)
@pytest.mark.parametrize("mu_b", [(3.3, 0.7)])
@pytest.mark.parametrize("hybridize", [True, False])
def test_laplace(mu_b: Tuple[float, float], hybridize: bool) -> None:
    '''
    Test to check that maximizing the likelihood recovers the parameters
    '''
    # test instance
    mu, b = mu_b

    # generate samples
    mus = mx.nd.zeros((NUM_SAMPLES, 1)) + mu
    bs = mx.nd.zeros((NUM_SAMPLES, 1)) + b

    laplace_distr = Laplace(mu=mus, b=bs)
    samples = laplace_distr.sample()

    init_biases = [
        mu - START_TOL_MULTIPLE * TOL * mu,
        inv_softplus(b + START_TOL_MULTIPLE * TOL * b),
    ]

    mu_hat, b_hat = maximum_likelihood_estimate_sgd(
        LaplaceOutput(), samples, hybridize=hybridize, init_biases=init_biases
    )

    assert (
        np.abs(mu_hat - mu) < TOL * mu
    ), f"mu did not match: mu = {mu}, mu_hat = {mu_hat}"
    assert (
        np.abs(b_hat - b) < TOL * b
    ), f"b did not match: b = {b}, b_hat = {b_hat}"

@pytest.mark.parametrize(
    "gamma, slopes, knot_spacings",
    [(2.0, np.array([3, 1, 3, 4]), np.array([0.3, 0.2, 0.35, 0.15]))],
)
@pytest.mark.parametrize("hybridize", [True, False])
def test_piecewise_linear(
    gamma: float,
    slopes: np.ndarray,
    knot_spacings: np.ndarray,
    hybridize: bool,
) -> None:
    '''
    Test to check that minimizing the CRPS recovers the quantile function
    '''
    num_samples = 500  # use a few samples for timeout failure

    gammas = mx.nd.zeros((num_samples,)) + gamma
    slopess = mx.nd.zeros((num_samples, len(slopes))) + mx.nd.array(slopes)
    knot_spacingss = mx.nd.zeros(
        (num_samples, len(knot_spacings))
    ) + mx.nd.array(knot_spacings)

    pwl_sqf = PiecewiseLinear(gammas, slopess, knot_spacingss)

    samples = pwl_sqf.sample()

    # Parameter initialization
    gamma_init = gamma - START_TOL_MULTIPLE * TOL * gamma
    slopes_init = slopes - START_TOL_MULTIPLE * TOL * slopes
    knot_spacings_init = knot_spacings
    # We perturb knot spacings such that even after the perturbation they sum to 1.
    mid = len(slopes) // 2
    knot_spacings_init[:mid] = (
        knot_spacings[:mid] - START_TOL_MULTIPLE * TOL * knot_spacings[:mid]
    )
    knot_spacings_init[mid:] = (
        knot_spacings[mid:] + START_TOL_MULTIPLE * TOL * knot_spacings[mid:]
    )

    init_biases = [gamma_init, slopes_init, knot_spacings_init]

    # check if it returns original parameters of mapped
    gamma_hat, slopes_hat, knot_spacings_hat = maximum_likelihood_estimate_sgd(
        PiecewiseLinearOutput(len(slopes)),
        samples,
        init_biases=init_biases,
        hybridize=hybridize,
        learning_rate=PositiveFloat(0.01),
        num_epochs=PositiveInt(20),
    )

    # Since the problem is highly non-convex we may not be able to recover the exact parameters
    # Here we check if the estimated parameters yield similar function evaluations at different quantile levels.
    quantile_levels = np.arange(0.1, 1.0, 0.1)

    # create a LinearSplines instance with the estimated parameters to have access to .quantile
    pwl_sqf_hat = PiecewiseLinear(
        mx.nd.array(gamma_hat),
        mx.nd.array(slopes_hat).expand_dims(axis=0),
        mx.nd.array(knot_spacings_hat).expand_dims(axis=0),
    )

    # Compute quantiles with the estimated parameters
    quantiles_hat = np.squeeze(
        pwl_sqf_hat.quantile(
            mx.nd.array(quantile_levels).expand_dims(axis=0), axis=1
        ).asnumpy()
    )

    # Compute quantiles with the original parameters
    # Since params is replicated across samples we take only the first entry
    quantiles = np.squeeze(
        pwl_sqf.quantile(
            mx.nd.array(quantile_levels)
            .expand_dims(axis=0)
            .repeat(axis=0, repeats=num_samples),
            axis=1,
        ).asnumpy()[0, :]
    )

    for ix, (quantile, quantile_hat) in enumerate(
        zip(quantiles, quantiles_hat)
    ):
        assert np.abs(quantile_hat - quantile) < TOL * quantile, (
            f"quantile level {quantile_levels[ix]} didn't match:"
            f" "
            f"q = {quantile}, q_hat = {quantile_hat}"
        )

@pytest.mark.skip("this test fails when run locally")
@pytest.mark.timeout(20)
@pytest.mark.parametrize("lambdas", [(0.1, 0.01)])
@pytest.mark.parametrize("mu_sigma", [(-1.5, 0.5)])
@pytest.mark.parametrize("hybridize", [True])
def test_box_cox_tranform(
    lambdas: Tuple[float, float],
    mu_sigma: Tuple[float, float],
    hybridize: bool,
):
    '''
    Test to check that maximizing the likelihood recovers the parameters
    '''
    # test instance
    lam_1, lam_2 = lambdas
    mu, sigma = mu_sigma

    # generate samples
    lamdas_1 = mx.nd.zeros((NUM_SAMPLES,)) + lam_1
    lamdas_2 = mx.nd.zeros((NUM_SAMPLES,)) + lam_2
    transform = InverseBoxCoxTransform(lamdas_1, lamdas_2)

    mus = mx.nd.zeros((NUM_SAMPLES,)) + mu
    sigmas = mx.nd.zeros((NUM_SAMPLES,)) + sigma
    gausian_distr = Gaussian(mus, sigmas)

    # Here the base distribution is Guassian which is transformed to
    # non-Gaussian via the inverse Box-Cox transform.
    # Sampling from `trans_distr` gives non-Gaussian samples
    trans_distr = TransformedDistribution(gausian_distr, transform)

    # Given the non-Gaussian samples find the true parameters
    # of the Box-Cox transformation as well as the underlying Gaussian distribution.
    samples = trans_distr.sample()

    init_biases = [
        mu - START_TOL_MULTIPLE * TOL * mu,
        inv_softplus(sigma - START_TOL_MULTIPLE * TOL * sigma),
        lam_1 - START_TOL_MULTIPLE * TOL * lam_1,
        inv_softplus(lam_2 - START_TOL_MULTIPLE * TOL * lam_2),
    ]

    mu_hat, sigma_hat, lam_1_hat, lam_2_hat = maximum_likelihood_estimate_sgd(
        TransformedDistributionOutput(
            GaussianOutput(),
            InverseBoxCoxTransformOutput(lb_obs=lam_2, fix_lambda_2=True),
        ),
        samples,
        init_biases=init_biases,
        hybridize=hybridize,
        learning_rate=PositiveFloat(0.01),
        num_epochs=PositiveInt(18),
    )

    assert (
        np.abs(lam_1_hat - lam_1) < TOL * lam_1
    ), f"lam_1 did not match: lam_1 = {lam_1}, lam_1_hat = {lam_1_hat}"
    # assert (
    #     np.abs(lam_2_hat - lam_2) < TOL * lam_2
    # ), f"lam_2 did not match: lam_2 = {lam_2}, lam_2_hat = {lam_2_hat}"

    assert np.abs(mu_hat - mu) < TOL * np.abs(
        mu
    ), f"mu did not match: mu = {mu}, mu_hat = {mu_hat}"
    assert (
        np.abs(sigma_hat - sigma) < TOL * sigma
    ), f"sigma did not match: sigma = {sigma}, sigma_hat = {sigma_hat}"

@pytest.mark.parametrize("num_bins", [6])
@pytest.mark.parametrize(
    "bin_probabilites", [np.array([0.3, 0.1, 0.05, 0.2, 0.1, 0.25])]
)
@pytest.mark.parametrize(
    "hybridize", [False]
)  # TODO: hybridize does not work in the test
def test_binned_likelihood(
    num_bins: float, bin_probabilites: np.ndarray, hybridize: bool
):
    '''
    Test to check that maximizing the likelihood recovers the parameters
    '''

    bin_prob = mx.nd.array(bin_probabilites)
    bin_center = mx.nd.array(np.logspace(-1, 1, num_bins))

    # generate samples
    bin_probs = mx.nd.zeros((NUM_SAMPLES, num_bins)) + bin_prob
    bin_centers = mx.nd.zeros((NUM_SAMPLES, num_bins)) + bin_center

    distr = Binned(bin_probs, bin_centers)
    samples = distr.sample()

    # add some jitter to the uniform initialization and normalize
    bin_prob_init = mx.nd.random_uniform(1 - TOL, 1 + TOL, num_bins) * bin_prob
    bin_prob_init = bin_prob_init / bin_prob_init.sum()

    init_biases = [bin_prob_init]

    bin_prob_hat, _ = maximum_likelihood_estimate_sgd(
        BinnedOutput(list(bin_center.asnumpy())),
        samples,
        init_biases=init_biases,
        hybridize=hybridize,
        learning_rate=PositiveFloat(0.05),
        num_epochs=PositiveInt(25),
    )

    assert all(
        mx.nd.abs(mx.nd.array(bin_prob_hat) - bin_prob) < TOL * bin_prob
    ), f"bin_prob did not match: bin_prob = {bin_prob}, bin_prob_hat = {bin_prob_hat}"

max_stars_repo_name: RubensZimbres/gluon-ts

max_issues_repo_name: RubensZimbres/gluon-ts

max_forks_repo_name: RubensZimbres/gluon-ts

reason: import numpy


text:
## The functions used to compare surface salinity
from __future__ import division, print_function
from salishsea_tools import (nc_tools,viz_tools,stormtools,tidetools)
from salishsea_tools.nowcast import figures
from datetime import datetime, timedelta
from pylab import *
#from sklearn import linear_model
from glob import glob
from IPython.core.display import HTML
from salishsea_tools.nowcast import figures
import matplotlib.pyplot as plt
import scipy.io as sio
import netCDF4 as nc
import numpy as np
import seaborn as sns
import math
#import glob
import os
import datetime

sns.set_style('darkgrid')

title_font = {
    'fontname': 'Bitstream Vera Sans', 'size': '15', 'color': 'black',
    'weight': 'medium'
}
axis_font = {'fontname': 'Bitstream Vera Sans', 'size': '13'}

ferry_stations = {'Tsawwassen': {'lat': 49.0084,'lon': -123.1281},
                  'Duke': {'lat': 49.1632,'lon': -123.8909},
                  'Vancouver': {'lat': 49.2827,'lon': -123.1207}}
paths = {'nowcast': '/data/dlatorne/MEOPAR/SalishSea/nowcast/',
        'longerresult': '/ocean/jieliu/research/meopar/river-treatment/14days_norefraserxml/',
        'widenresult': '/data/jieliu/MEOPAR/river-treatment/24nor_NW/' }

def results_dataset(period, grid, results_dir):
    """Return the results dataset for period (e.g. 1h or 1d)
    and grid (e.g. grid_T, grid_U) from results_dir.
    """
    filename_pattern = 'SalishSea_{period}_*_{grid}.nc'
    filepaths = glob(os.path.join(results_dir, filename_pattern.format(period=period, grid=grid)))
    return  nc.Dataset(filepaths[0])

def date(year, month, day, day_start, day_end, results_home,  period, grid):

    day_range = np.arange(day_start, day_end+1)
    day_len = len(day_range)
    files_all = [None] * day_len
    inds = np.arange(day_len)

    for i, day in zip(inds, day_range):
        run_date = datetime.datetime(year,month, day)
        #results_home = '/data/jieliu/MEOPAR/river-treatment/24nor_NW/'
        results_dir = os.path.join(results_home, run_date.strftime('%d%b%y').lower())
        filename = 'SalishSea_' + period + '_' + run_date.strftime('%Y%m%d').lower() + \
        '_' + run_date.strftime('%Y%m%d').lower() + '_' + grid + '.nc'
        file_single = os.path.join(results_dir, filename)
        files_all[i] = file_single

    return files_all

def find_dist (q, lon11, lat11, X, Y, bathy, longitude, latitude, saline_nemo_3rd, saline_nemo_4rd):
    k=0
    values =0
    valuess=0
    dist = np.zeros(9)
    weights = np.zeros(9)
    value_3rd=np.zeros(9)
    value_4rd=np.zeros(9)
    #regr =linear_model.LinearRegression()
    #regr.fit(lon11,lat11);
    #regr.coef_

    [x1, j1] = tidetools.find_closest_model_point(lon11[q],lat11[q]),\
                                        X,Y,bathy,lon_tol=0.0052,lat_tol=0.00210,allow_land=False)
    for i in np.arange(x1-1,x1+2):
        for j in np.arange(j1-1,j1+2):
            dist[k]=tidetools.haversine(lon11[q],lat11[q],longitude[i,j],latitude[i,j])
            weights[k]=1.0/dist[k]
            value_3rd[k]=saline_nemo_3rd[i,j]*weights[k]
            value_4rd[k]=saline_nemo_4rd[i,j]*weights[k]
            values=values+value_3rd[k]
            valuess=valuess+value_4rd[k]
            k+=1

    return values, valuess, weights

def find_dist_ave (q, lon11, lat11, X, Y, bathy, longitude, latitude, saline_nemo_3rd, saline_nemo_4rd):
    k=0
    values_ave =0
    valuess_ave=0
    dist = np.zeros(9)
    weights = np.zeros(9)
    value_3rd=np.zeros(9)
    value_4rd=np.zeros(9)
    #regr =linear_model.LinearRegression()
    #regr.fit(lon11,lat11);
    #regr.coef_

    [x1, j1] = tidetools.find_closest_model_point(lon11[q],lat11[q]),\
                                        X,Y,bathy,lon_tol=0.0052,lat_tol=0.00210,allow_land=False)
    for i in np.arange(x1-1,x1+2):
        for j in np.arange(j1-1,j1+2):
            dist[k]=tidetools.haversine(lon11[q],lat11[q],longitude[i,j],latitude[i,j])
            weights[k]=1.0/dist[k]
            value_3rd[k]=saline_nemo_3rd[i,j]*weights[k]
            value_4rd[k]=saline_nemo_4rd[i,j]*weights[k]
            values_ave=values_ave+value_3rd[k]
            valuess_ave=valuess_ave+value_4rd[k]
            k+=1

    return values_ave, valuess_ave, weights

def get_SS2_bathy_data():
    """Get the original Salish Sea 2 bathymetry and grid data
    e.g. bathy, X, Y = get_SS2_bathy_data()

    .. note::

        This function is deprecated due to hard-coding of
        :file:`/ocean/klesouef/` path.
        Use :py:func:`tidetools.get_bathy_data` instead.

    :returns: bathy, X, Y
    """
    grid = nc.Dataset(
        '/ocean/jieliu/research/meopar/nemo-forcing/grid/bathy_meter_SalishSea2.nc', 'r')
    bathy = grid.variables['Bathymetry'][:, :]
    X = grid.variables['nav_lon'][:, :]
    Y = grid.variables['nav_lat'][:, :]
    return bathy, X, Y

def get_SS5_bathy_data():
    """Get the Salish Sea 5 bathymetry and grid data
    e.g. bathy, X, Y = get_SS5_bathy_data()

    .. note::

        This function is deprecated due to hard-coding of
        :file:`/ocean/klesouef/` path.
        Use :py:func:`tidetools.get_bathy_data` instead.

    :returns: bathy, X, Y
    """
    grid = nc.Dataset(
        '/ocean/jieliu/research/meopar/river-treatment/bathy_meter_SalishSea5.nc', 'r')
    bathy = grid.variables['Bathymetry'][:, :]
    X = grid.variables['nav_lon'][:, :]
    Y = grid.variables['nav_lat'][:, :]
    return bathy, X, Y

def get_SS6_bathy_data():
    """Get the Salish Sea 6 bathymetry and grid data
    e.g. bathy, X, Y = get_SS6_bathy_data()

    .. note::

        This function is deprecated due to hard-coding of
        :file:`/ocean/klesouef/` path.
        Use :py:func:`tidetools.get_bathy_data` instead.

    :returns: bathy, X, Y
    """
    grid = nc.Dataset(
        '/ocean/jieliu/research/meopar/river-treatment/bathy_meter_SalishSea6.nc', 'r')
    bathy = grid.variables['Bathymetry'][:, :]
    X = grid.variables['nav_lon'][:, :]
    Y = grid.variables['nav_lat'][:, :]
    return bathy, X, Y

def salinity_fxn(saline, run_date, results_home):
    """The significance of this function was to return longitude,
       latitude, salinity values for observations, 1.5m of 3rd & 4rd,
       3m average of 3rd & 4rd model result """

    a=saline['ferryData']
    b=a['data']
    dataa = b[0,0]
    time=dataa['matlabtime'][0,0]
    lonn=dataa['Longitude'][0,0]
    latt=dataa['Latitude'][0,0]
    salinity=dataa['Practical_Salinity'][0,0]

    a=len(time)
    lon1=np.zeros([a,1])
    lat1=np.zeros([a,1])
    salinity1=np.zeros([a,1])
    for i in np.arange(0,a):
        matlab_datenum = np.float(time[i])
        python_datetime = datetime.datetime.fromordinal(int(matlab_datenum))\
        + timedelta(days=matlab_datenum%1) - timedelta(days = 366)

        if((python_datetime.year == run_date.year) & (python_datetime.month == run_date.month)\
           & (python_datetime.day == run_date.day)
           & (python_datetime.hour >= 3))&(python_datetime.hour < 5):
            lon1[i]=lonn[i]
            lat1[i]=latt[i]
            salinity1[i]=salinity[i]

    mask=lon1[:,0]!=0
    lon1_2_4=lon1[mask]
    lat1_2_4=lat1[mask]
    salinity1_2_4=salinity1[mask]
    lon11=lon1_2_4[0:-1:20]
    lat11=lat1_2_4[0:-1:20]
    salinity11=salinity1_2_4[0:-1:20]
    if results_home == paths['longerresult']:
        bathynew, X, Y = get_SS5_bathy_data()
    elif results_home == paths['nowcast']:
        bathyold, X, Y = get_SS2_bathy_data()

    #bathy, X, Y = tidetools.get_SS2_bathy_data()

    #aa=date(run_date.year,run_date.month,run_date.day,run_date.day,'1h','grid_T')
    #sim_date = datetime.datetime(2015,3,19)####need to change for \
    #different daily model results, construct a datetime object
    #run_date = datetime.datetime(2015,3,19)

    date_str = run_date.strftime('%d-%b-%Y') ##create a string based on this date
    filepath_name = date(run_date.year,run_date.month, run_date.day,\
    run_date.day,run_date.day, results_home,'1h','grid_T')
    tracers=nc.Dataset(filepath_name[0])
    #j=int(aa[0][65:67])
    #jj=int(aa[0][67:69])
    latitude=tracers.variables['nav_lat'][:]
    longitude=tracers.variables['nav_lon'][:]
    saline_nemo = tracers.variables['vosaline']
    saline_nemo_3rd = saline_nemo[3,1, 0:898, 0:398]
    saline_nemo_4rd = saline_nemo[4,1, 0:898, 0:398]
    saline_nemo_ave3rd = np.mean(saline_nemo[3, 0:3, 0:898, 0:398], axis = 0)
    saline_nemo_ave4rd = np.mean(saline_nemo[4, 0:3, 0:898, 0:398], axis = 0)

    matrix=np.zeros([36,9])
    matrix_ave=np.zeros([36,9])
    values=np.zeros([36,1])
    valuess=np.zeros([36,1])
    values_ave=np.zeros([36,1])
    valuess_ave=np.zeros([36,1])
    value_mean_3rd_hour=np.zeros([36,1])
    value_mean_ave3rd=np.zeros([36,1])
    value_mean_ave4rd=np.zeros([36,1])
    value_mean_4rd_hour=np.zeros([36,1])
    for q in np.arange(0,36):
        if results_home == paths['longerresult']:
            values[q], valuess[q], matrix[q,:]=find_dist(q, lon11, lat11, X, Y,\
                                     bathynew, longitude, latitude, saline_nemo_3rd, saline_nemo_4rd)
            value_mean_3rd_hour[q]=values[q]/sum(matrix[q])
            value_mean_4rd_hour[q]=valuess[q]/sum(matrix[q])

            values_ave[q], valuess_ave[q], matrix_ave[q,:]=find_dist_ave(q, lon11, lat11, X, Y,\
                                     bathynew, longitude, latitude, saline_nemo_ave3rd, saline_nemo_ave4rd)
            value_mean_ave3rd[q]=values_ave[q]/sum(matrix_ave[q])
            value_mean_ave4rd[q]=valuess_ave[q]/sum(matrix_ave[q])
        elif results_home == paths['nowcast']:
            values[q], valuess[q], matrix[q,:]=find_dist(q, lon11, lat11, X, Y,\
                                     bathyold, longitude, latitude, saline_nemo_3rd, saline_nemo_4rd)
            value_mean_3rd_hour[q]=values[q]/sum(matrix[q])
            value_mean_4rd_hour[q]=valuess[q]/sum(matrix[q])

            values_ave[q], valuess_ave[q], matrix_ave[q,:]=find_dist_ave(q, lon11, lat11, X, Y,\
                                     bathyold, longitude, latitude, saline_nemo_ave3rd, saline_nemo_ave4rd)
            value_mean_ave3rd[q]=values_ave[q]/sum(matrix_ave[q])
            value_mean_ave4rd[q]=valuess_ave[q]/sum(matrix_ave[q])

    return lon11, lat11, lon1_2_4, lat1_2_4,\
    value_mean_3rd_hour, value_mean_4rd_hour, \
    value_mean_ave3rd, value_mean_ave4rd,\
    salinity11, salinity1_2_4, date_str

def salinity_fxn_five_times(saline, run_date, results_home):
    """The significance of this function was to return longitude,
       latitude, salinity values for observations, 1.5m of 3rd & 4rd,
       3m average of 3rd & 4rd model result """

    a=saline['ferryData']
    b=a['data']
    dataa = b[0,0]
    time=dataa['matlabtime'][0,0]
    lonn=dataa['Longitude'][0,0]
    latt=dataa['Latitude'][0,0]
    salinity=dataa['Practical_Salinity'][0,0]

    a=len(time)
    lon1=np.zeros([a,1])
    lat1=np.zeros([a,1])
    salinity1=np.zeros([a,1])
    for i in np.arange(0,a):
        matlab_datenum = np.float(time[i])
        python_datetime = datetime.datetime.fromordinal(int(matlab_datenum))\
        + timedelta(days=matlab_datenum%1) - timedelta(days = 366)

        if((python_datetime.year == run_date.year) & (python_datetime.month == run_date.month)\
           & (python_datetime.day == run_date.day)
           & (python_datetime.hour >= 3))&(python_datetime.hour < 5):
            lon1[i]=lonn[i]
            lat1[i]=latt[i]
            salinity1[i]=salinity[i]

    mask=lon1[:,0]!=0
    lon1_2_4=lon1[mask]
    lat1_2_4=lat1[mask]
    salinity1_2_4=salinity1[mask]
    lon11=lon1_2_4[0:-1:20]
    lat11=lat1_2_4[0:-1:20]
    salinity11=salinity1_2_4[0:-1:20]
    if results_home == paths['longerresult']:
        bathynew, X, Y = get_SS5_bathy_data()
    elif results_home == paths['nowcast']:
        bathyold, X, Y = get_SS2_bathy_data()

    date_str = run_date.strftime('%d-%b-%Y') ##create a string based on this date
    filepath_name = date(run_date.year,run_date.month, run_date.day,\
    run_date.day,run_date.day, results_home,'1h','grid_T')
    tracers=nc.Dataset(filepath_name[0])
    latitude=tracers.variables['nav_lat'][:]
    longitude=tracers.variables['nav_lon'][:]
    saline_nemo = tracers.variables['vosaline']
    saline_nemo_3rd = saline_nemo[3,1, 0:898, 0:398]
    saline_nemo_4rd = saline_nemo[4,1, 0:898, 0:398]
    saline_nemo_ave3rd = np.mean(saline_nemo[3, 0:3, 0:898, 0:398], axis = 0)
    saline_nemo_ave4rd = np.mean(saline_nemo[4, 0:3, 0:898, 0:398], axis = 0)

    matrix=np.zeros([len(lon11),9])
    matrix_ave=np.zeros([len(lon11),9])
    values=np.zeros([len(lon11),1])
    valuess=np.zeros([len(lon11),1])
    values_ave=np.zeros([len(lon11),1])
    valuess_ave=np.zeros([len(lon11),1])
    value_mean_3rd_hour=np.zeros([len(lon11),1])
    value_mean_ave3rd=np.zeros([len(lon11),1])
    value_mean_ave4rd=np.zeros([len(lon11),1])
    value_mean_4rd_hour=np.zeros([len(lon11),1])
    for q in np.arange(0,len(lon11)):
        if results_home == paths['longerresult']:
            values[q], valuess[q], matrix[q,:]=find_dist(q, lon11, lat11, X, Y,\
                                     bathynew, longitude, latitude, saline_nemo_3rd, saline_nemo_4rd)
            value_mean_3rd_hour[q]=values[q]/sum(matrix[q])
            value_mean_4rd_hour[q]=valuess[q]/sum(matrix[q])

            values_ave[q], valuess_ave[q], matrix_ave[q,:]=find_dist_ave(q, lon11, lat11, X, Y,\
                                     bathynew, longitude, latitude, saline_nemo_ave3rd, saline_nemo_ave4rd)
            value_mean_ave3rd[q]=values_ave[q]/sum(matrix_ave[q])
            value_mean_ave4rd[q]=valuess_ave[q]/sum(matrix_ave[q])
        elif results_home == paths['nowcast']:
            values[q], valuess[q], matrix[q,:]=find_dist(q, lon11, lat11, X, Y,\
                                     bathyold, longitude, latitude, saline_nemo_3rd, saline_nemo_4rd)
            value_mean_3rd_hour[q]=values[q]/sum(matrix[q])
            value_mean_4rd_hour[q]=valuess[q]/sum(matrix[q])

            values_ave[q], valuess_ave[q], matrix_ave[q,:]=find_dist_ave(q, lon11, lat11, X, Y,\
                                     bathyold, longitude, latitude, saline_nemo_ave3rd, saline_nemo_ave4rd)
            value_mean_ave3rd[q]=values_ave[q]/sum(matrix_ave[q])
            value_mean_ave4rd[q]=valuess_ave[q]/sum(matrix_ave[q])

    return lon11, lat11, lon1_2_4, lat1_2_4,\
    value_mean_3rd_hour, value_mean_4rd_hour, \
    value_mean_ave3rd, value_mean_ave4rd,\
    salinity11, salinity1_2_4, date_str

#def j_i_model_points(saline, run_date, results_home):
    #"""This function was made to find the ferry route longitude and
	#latitude for model results, return lon and lat indices in list"""
    #lon11, lat11, lon1_2_4, lat1_2_4,\
    #value_mean_3rd_hour, value_mean_4rd_hour, \
    #value_mean_ave3rd, value_mean_ave4rd,\
    #salinity11, salinity1_2_4,date_str = salinity_fxn(saline, run_date, results_home)

    #[x1, j1] = tidetools.find_closest_model_point(lon11[q],lat11[q]),\
                                        #X,Y,bathy,lon_tol=0.0052,lat_tol=0.00210,allow_land=False)

def salinity_ferry_route(grid_T, grid_B, PNW_coastline, sal_hr, saline, run_date,results_home):
    """ plot daily salinity comparisons between ferry observations
    and model results as well as ferry route with model salinity
    distribution.

    :arg grid_B: Bathymetry dataset for the Salish Sea NEMO model.
    :type grid_B: :class:`netCDF4.Dataset`

    :arg PNW_coastline: Coastline dataset.
    :type PNW_coastline: :class:`mat.Dataset`

    :arg ferry_sal: saline
    :type ferry_sal: numpy

    :returns: fig
    """

    latitude=grid_T.variables['nav_lat']
    longitude=grid_T.variables['nav_lon']
    fig, axs = plt.subplots(1, 2, figsize=(15, 8))

    figures.plot_map(axs[1], grid_B, PNW_coastline)
    axs[1].set_xlim(-124.5, -122.5)
    axs[1].set_ylim(48.2, 49.5)
    viz_tools.set_aspect(axs[1],coords='map',lats=latitude)
    cmap=plt.get_cmap('spectral')
    cmap.set_bad('burlywood')
    mesh=axs[1].pcolormesh(longitude[:],latitude[:],sal_hr[:],cmap=cmap)
    cbar=fig.colorbar(mesh)
    plt.setp(plt.getp(cbar.ax.axes, 'yticklabels'), color='w')
    cbar.set_label('Pratical Salinity', color='white')

    axs[1].set_title('Ferry Route: 3am[UTC] 1.5m model result ', **title_font)

    bbox_args = dict(boxstyle='square', facecolor='white', alpha=0.7)
    stations=['Tsawwassen','Duke','Vancouver']
    for stn in stations:
        axs[1].plot(ferry_stations[stn]['lon'], ferry_stations[stn]['lat'], marker='D', \
                    color='white',\
                 markersize=10, markeredgewidth=2)
    axs[1].annotate ('Tsawwassen',(ferry_stations['Tsawwassen']['lon'] + 0.02,\
    ferry_stations['Tsawwassen']['lat'] + 0.12), fontsize=15, color='black', bbox=bbox_args )
    axs[1].annotate ('Duke',(ferry_stations['Duke']['lon'] - 0.35,\
    ferry_stations['Duke']['lat'] ),fontsize=15, color='black', bbox=bbox_args )
    axs[1].annotate ('Vancouver',(ferry_stations['Vancouver']['lon'] - 0.1,\
    ferry_stations['Vancouver']['lat']+ 0.09 ),fontsize=15, color='black', bbox=bbox_args )
    figures.axis_colors(axs[1], 'white')

    #filepath_name = date(run_date.year,run_date.month, run_date.day,\
     #run_date.day,run_date.day, results_home,'1h','grid_T')
    lon11, lat11, lon1_2_4, lat1_2_4,\
    value_mean_3rd_hour, value_mean_4rd_hour, \
    value_mean_ave3rd, value_mean_ave4rd,\
    salinity11, salinity1_2_4,date_str = salinity_fxn(saline, run_date, results_home )
    axs[1].plot(lon11,lat11,'black', linewidth = 4)
    model_salinity_3rd_hour=axs[0].plot(lon11,value_mean_3rd_hour,'DodgerBlue',\
                                    linewidth=2, label='3 am [UTC]')
    model_salinity_4rd_hour=axs[0].plot(lon11,value_mean_4rd_hour,'MediumBlue',\
                                        linewidth=2, label="4 am [UTC]" )
    observation_salinity=axs[0].plot(lon1_2_4,salinity1_2_4,'DarkGreen', \
                                     linewidth=2, label="Observed")
    axs[0].text(0.25, -0.1,'Observations from Ocean Networks Canada', \
                transform=axs[0].transAxes, color='white')

    axs[0].set_xlim(-124, -123)
    axs[0].set_ylim(0, 30)
    axs[0].set_title('Surface Salinity: ' + date_str, **title_font)
    axs[0].set_xlabel('Longitude', **axis_font)
    axs[0].set_ylabel('Practical Salinity', **axis_font)
    axs[0].legend()
    axs[0].grid()

    fig.patch.set_facecolor('#2B3E50')
    figures.axis_colors(axs[0], 'gray')

    return fig

def salinity_ferry_route_more(grid_T, grid_B, PNW_coastline, ave, sal_hr_1, sal_hr_ave, saline, run_date,results_home):
    """ plot daily salinity comparisons between ferry observations
    and model results as well as ferry route with model salinity
    distribution.

    :arg grid_B: Bathymetry dataset for the Salish Sea NEMO model.
    :type grid_B: :class:`netCDF4.Dataset`

    :arg PNW_coastline: Coastline dataset.
    :type PNW_coastline: :class:`mat.Dataset`

    :arg ferry_sal: saline
    :type ferry_sal: numpy

    :returns: fig
    """

    latitude=grid_T.variables['nav_lat']
    longitude=grid_T.variables['nav_lon']
    fig, axs = plt.subplots(1, 2, figsize=(15, 8))

    figures.plot_map(axs[1], grid_B, PNW_coastline)
    axs[1].set_xlim(-124.5, -122.5)
    axs[1].set_ylim(48.2, 49.5)
    viz_tools.set_aspect(axs[0],coords='map',lats=latitude)
    cmap=plt.get_cmap('spectral')
    cmap.set_bad('burlywood')
    if ave ==0:
        mesh=axs[1].pcolormesh(longitude[:],latitude[:],sal_hr_1[:],cmap=cmap)
        axs[1].set_title('Ferry Route: 3am[UTC] 1.5m model result ', **title_font)
    elif ave ==1:
        mesh=axs[1].pcolormesh(longitude[:],latitude[:],sal_hr_ave[:],cmap=cmap)
        axs[1].set_title('Ferry Route: 3am[UTC] 3m average model result ', **title_font)
    cbar=fig.colorbar(mesh)
    plt.setp(plt.getp(cbar.ax.axes, 'yticklabels'), color='w')
    cbar.set_label('Pratical Salinity', color='white')

    bbox_args = dict(boxstyle='square', facecolor='white', alpha=0.7)
    stations=['Tsawwassen','Duke','Vancouver']
    for stn in stations:
        axs[1].plot(ferry_stations[stn]['lon'], ferry_stations[stn]['lat'], marker='D', \
                    color='white',\
                 markersize=10, markeredgewidth=2)
    axs[1].annotate ('Tsawwassen',(ferry_stations['Tsawwassen']['lon'] + 0.02,\
    ferry_stations['Tsawwassen']['lat'] + 0.12), fontsize=15, color='black', bbox=bbox_args )
    axs[1].annotate ('Duke',(ferry_stations['Duke']['lon'] - 0.35,\
    ferry_stations['Duke']['lat'] ),fontsize=15, color='black', bbox=bbox_args )
    axs[1].annotate ('Vancouver',(ferry_stations['Vancouver']['lon'] - 0.1,\
    ferry_stations['Vancouver']['lat']+ 0.09 ),fontsize=15, color='black', bbox=bbox_args )
    figures.axis_colors(axs[1], 'white')

    #filepath_name = date(run_date.year,run_date.month, run_date.day,\
     #run_date.day,run_date.day, results_home,'1h','grid_T')
    lon11, lat11, lon1_2_4, lat1_2_4,\
    value_mean_3rd_hour, value_mean_4rd_hour, \
    value_mean_ave3rd, value_mean_ave4rd,\
    salinity11, salinity1_2_4,date_str = salinity_fxn(saline, run_date, results_home)
    axs[1].plot(lon11,lat11,'black', linewidth = 4)
    if ave ==0:
        model_salinity_3rd_hour=axs[0].plot(lon11,value_mean_3rd_hour,'DodgerBlue',\
                                    linewidth=2, label='3 am [UTC]')
        model_salinity_4rd_hour=axs[0].plot(lon11,value_mean_4rd_hour,'MediumBlue',\
                                        linewidth=2, label="4 am [UTC]" )
        observation_salinity=axs[0].plot(lon1_2_4,salinity1_2_4,'DarkGreen', \
                                     linewidth=2, label="Observed")
    elif ave ==1:
        model_salinity_3rd_hour=axs[0].plot(lon11,value_mean_ave3rd,'DodgerBlue',\
                                    linewidth=2, label='3 am [UTC]')
        model_salinity_4rd_hour=axs[0].plot(lon11,value_mean_ave4rd,'MediumBlue',\
                                        linewidth=2, label="4 am [UTC]" )
        observation_salinity=axs[0].plot(lon1_2_4,salinity1_2_4,'DarkGreen', \
                                     linewidth=2, label="Observed")
    axs[0].text(0.25, -0.1,'Observations from Ocean Networks Canada', \
                transform=axs[0].transAxes, color='white')

    axs[0].set_xlim(-124, -123)
    axs[0].set_ylim(0, 30)
    axs[0].set_title('Surface Salinity: ' + date_str, **title_font)
    axs[0].set_xlabel('Longitude', **axis_font)
    axs[0].set_ylabel('Practical Salinity', **axis_font)
    axs[0].legend()
    axs[0].grid()
    fig.patch.set_facecolor('#2B3E50')
    figures.axis_colors(axs[0], 'gray')

    return fig

def find_min_value_location(run_date, results_home, saline):
    """This function was made to find out the minimum salinity value
       and longitude for observation and nowcasts(1.5m and 3m ave) or
       for observation and new results(1.5m and 3m ave)."""

    #filepath_name = date(run_date.year,run_date.month, run_date.day,\
     #run_date.day,run_date.day, results_home,'1h','grid_T')
    lon11, lat11, lon1_2_4, lat1_2_4,\
    value_mean_3rd_hour, value_mean_4rd_hour, \
    value_mean_ave3rd, value_mean_ave4rd,\
    salinity11, salinity1_2_4,date_str = salinity_fxn(saline, run_date, results_home)

    ##For observation
    salinity_min_obs = nanmin(salinity1_2_4)
    ind_obs = nanargmin(salinity1_2_4)
    lon_min_obs = lon1_2_4[ind_obs]

    ## For nowcast or new model result with 1.5m depth for 3rd model time
    salinity_min_15 = nanmin(value_mean_3rd_hour)
    ind_15 = nanargmin(value_mean_3rd_hour)
    lon_min_15 = lon11[ind_15]

    ## For nowcast or new model result with average 3m depth for 3rd model time
    salinity_min_ave = nanmin(value_mean_ave3rd)
    ind_ave = nanargmin(value_mean_ave3rd)
    lon_min_ave = lon11[ind_ave]

    return salinity_min_obs, lon_min_obs, salinity_min_15,\
            lon_min_15, salinity_min_ave, lon_min_ave

def freshwater_amount(saline, run_date, results_home):
    """This function was to calculate freshwater amount for observation, 1.5m &
        average 3m model nowcasts & new results, we assume salinity below 30 is
        freshwater"""
    lon11, lat11, lon1_2_4, lat1_2_4,\
    value_mean_3rd_hour, value_mean_4rd_hour, \
    value_mean_ave3rd, value_mean_ave4rd,\
    salinity11, salinity1_2_4,date_str = salinity_fxn(saline, run_date, results_home)

    salinity1_2_4[np.isnan(salinity1_2_4)] = 30 ##set nan to 30 psu for observation

    obs_Sdx = np.zeros(len(lon1_2_4)-1) ## index from 0 to 718
    mod_now_15Sdx = np.zeros(len(lon11)-1)
    mod_new_15Sdx = np.zeros(len(lon11)-1)
    mod_now_aveSdx = np.zeros(len(lon11)-1)
    mod_new_aveSdx = np.zeros(len(lon11)-1)

    length_obs = np.arange(len(lon1_2_4)-1)
    length_mod = np.arange(len(lon11) - 1)
    ## For observation:
    for i in length_obs:
        obs_Sdx[i] = (30 - salinity1_2_4[:,0][i]) * -np.diff(lon1_2_4[:,0])[i] * \
        111000 * cos(pi*lat1_2_4[i]/180) #discrete outcome by multiplying
    obs_total_integral = np.cumsum(obs_Sdx)
    max_amount_obs = max(obs_total_integral)
    if max_amount_obs <= 0:
        max_amount_obs == -min(obs_total_integral)
    ## For nowcasts results:
    if results_home == paths['nowcast']:    #fresh water salinity(diff between 28 and salinity) and dx
        for j in length_mod:
            mod_now_15Sdx[j] = (30 - value_mean_3rd_hour[:,0][j]) * -np.diff(lon11[:,0])[j] * \
            111000 * cos(pi*lat11[j]/180)  ## whether np.diff is positive or not depends
                                      # on the ferry route
            mod_now_aveSdx[j] = (30 - value_mean_ave3rd[:,0][j]) * -np.diff(lon11[:,0])[j] * \
            111000 * cos(pi*lat11[j]/180)
        mod_total_now_15integral = np.cumsum(mod_now_15Sdx)
        mod_total_now_aveintegral = np.cumsum(mod_now_aveSdx)
        ## Find max for 1.5m & average 3m depth nowcast result
        max_mod_now15 = max(mod_total_now_15integral)
        max_mod_nowave = max(mod_total_now_aveintegral)

        if max_mod_now15 <= 0:
            max_mod_now15 == -min(mod_total_now_15integral)
        if max_mod_nowave <=0:
            max_mod_nowave == -min(mod_total_now_aveintegral)

        return max_amount_obs, max_mod_now15, max_mod_nowave

    if results_home == paths['longerresult']:
        for j in length_mod:
            mod_new_15Sdx[j] = (30 - value_mean_3rd_hour[:,0][j]) * -np.diff(lon11[:,0])[j] * \
            111000 * cos(pi*lat11[j]/180)
            mod_new_aveSdx[j] = (30 - value_mean_ave3rd[:,0][j]) * -np.diff(lon11[:,0])[j] * \
            111000 * cos(pi*lat11[j]/180)
        mod_total_new_15integral = np.cumsum(mod_new_15Sdx)
        mod_total_new_aveintegral = np.cumsum(mod_new_aveSdx)
        ## Find max for 1.5m & average 3m depth new result
        max_mod_new15 = max(mod_total_new_15integral)
        max_mod_newave = max(mod_total_new_aveintegral)

        if max_mod_new15 <=0:
            max_mod_new15 == -min(mod_total_new_15integral)
        if max_mod_newave <=0:
            max_mod_newave == -min(mod_total_new_aveintegral)
        return max_amount_obs, max_mod_new15, max_mod_newave

def plot_freshwater_amount(obs_amount, mod15_now_amount, mod15_new_amount,modave_now_amount, modave_new_amount):
    """This function was made to plot time series of amount of fresh water from June 16 to 29
        for observations and model"""
    fig, axs = plt.subplots(1, 2, figsize=(100, 40))
    ## for time defination
    time =[]
    for t in np.arange(11):
        time.append(t)

    ## xtick and xticklabels for plot
    group_labels = ['06/16','06/24','06/29']

    ## observation & 1.5m nowcast & new model result minimim values
    ax = axs[0]
    ax.plot(time, obs_amount,'b-', marker = 'o', markersize = 35,linewidth=5.0, label = 'observed value')
    ax.plot(time, mod15_now_amount,'g-',marker = '^', markersize = 35, linewidth=5.0, label = '1.5m nowcast value')
    ax.plot(time, mod15_new_amount,'y-', marker = 's', markersize = 35,linewidth=5.0, label = '1.5m new result value')
    plt.setp(ax, xticks=[0, 4, 10 ], xticklabels=group_labels)
    plt.setp(ax.get_xticklabels(), fontsize=65)
    plt.setp(ax.get_yticklabels(), fontsize=65)
    ax.set_title('Total freshwater amount of 1.5m depth ', fontsize = 80)
    #ax.set_xlim(0, 10)
    #ax.set_ylim(0, 20)
    ax.set_xlabel('Date', fontsize = 65)
    ax.set_ylabel('Total freshwater amount [m]', fontsize = 65)
    ax.grid('on')
    ax.legend(loc = 2, fontsize = 55)

    ## observation & average 3m nowcast & new model result minimim values
    ax = axs[1]
    ax.plot(time, obs_amount,'bo-', marker = 'o', markersize = 35,linewidth=5.0,label = 'observed value')
    ax.plot(time, modave_now_amount,'g-', marker = '^', markersize = 35,linewidth=5.0,label = 'average 3m nowcast value')
    ax.plot(time, modave_new_amount,'y-', marker = 's', markersize = 35,linewidth=5.0,label = 'average 3m new result value')
    plt.setp(ax, xticks=[0, 4, 10 ], xticklabels=group_labels)
    plt.setp(ax.get_xticklabels(), fontsize=65)
    plt.setp(ax.get_yticklabels(), fontsize=65)
    ax.set_title('Total freshwater amount of the average 3m depth ', fontsize = 80)
    ax.set_xlabel('Date', fontsize = 65)
    ax.set_ylabel('Total freshwater amount [m]', fontsize = 65)
    ax.grid('on')
    ax.legend(loc = 2,fontsize = 55)
    return fig

max_stars_repo_name: SalishSeaCast/analysis

max_issues_repo_name: SalishSeaCast/analysis

max_forks_repo_name: SalishSeaCast/analysis


text:
"""
File: gov-bills-scraper.py
Author: Zachary King
Description: A web scraper that retrieves the
	stock value for a given list of companies
	(by symbols), and plots the values on a
	bar graph accordingly, with color-coding.
"""

import requests
from bs4 import BeautifulSoup
import numpy as np
import matplotlib.pyplot as plt

# Variable declarations
symbols = ["aapl", "goog", "msft", "yhoo", "amzn", "fb", "twtr", "nflx"]
numCompanies = len(symbols)
values = []
arrows = []

# Scrape the data
for s in symbols:
	url = "http://finance.yahoo.com/q?&s=" + s
	r = requests.get(url)
	soup = BeautifulSoup(r.content, "html.parser")
	price = soup.find_all("span", {"id": "yfs_l84_"+s})
	arrow = soup.find_all("span", {"id": "yfs_c63_"+s})
	arrows.append(arrow[0].next_element["class"])
	print(s.upper() + " share value: $" + price[0].text)
	values.append(price[0].text)
	symbols[symbols.index(s)] = symbols[symbols.index(s)].upper()

# convert values from strings to floats
for i in range(numCompanies):
	values[i] = float(values[i])

# Plot and figure
fig = plt.figure()
fig.canvas.set_window_title('Stock-Scraper')
ax = fig.add_subplot(111)

index = np.arange(numCompanies)
rects = ax.bar(index, values, 0.25, alpha=0.5)

# Show the number value for the share above the bar
for i in range(numCompanies):
	ax.annotate('$' + str(values[i]), xy=(i+0.12, values[i]), xytext=(i - 0.12, values[i] + 10),
            arrowprops=dict(arrowstyle='-'),
            )

# Labels and titles on graph
plt.xlabel('Company')
plt.ylabel('Share Value')
plt.title('Current Share Value for Companies')
plt.xticks(index + 0.25, symbols)

# Color the bars based on the companies' stock status
for i, arrow in enumerate(arrows):
	if (arrow == [u'pos_arrow']):
		rects[i].set_color('g') # stock is going up
	elif (arrow == [u'neg_arrow']):
		rects[i].set_color('r') # stock is going down
	else:
		rects[i].set_color('b') # unknown data

# Show the graph
fig.tight_layout()
plt.show()

max_stars_repo_name: zach-king/CoolPython

max_issues_repo_name: zach-king/CoolPython

max_forks_repo_name: zach-king/CoolPython


text:
# ToDO:

import logging
import pickle
import string
import sys
import time

import numpy as np

import utils
from templates.template import TemplateBaseClass
from tqdm import tqdm
import multiprocessing
from IPython.core.debugger import Pdb
import os
#os.system("taskset -p 0xfffff %d" % os.getpid())

class Template6(TemplateBaseClass):
    """
    (e1 r e2) since (e1 r1 u1) and (u1 r2 e2) and ((r1.r2)~r)
    """

    def __init__(self, kblist, base_model, use_hard_triple_scoring=True, load_table=None, dump_file=None, parts=1, offset=0):
        super().__init__()
        self.kb = kblist[0]
        self.kb_val = kblist[1]
        self.kb_test = kblist[2]
        self.base_model = base_model
        self.use_hard_triple_scoring = use_hard_triple_scoring
        # Pdb().set_trace()
        #self.exp_template = '$html_fact_rprime'
        self.exp_template = '$html_fact_e1_r1_u1 and $html_fact_u1_r2_e2'
        self.parts = parts
        self.offset = offset
        self.dump_file = dump_file
        if(load_table == None):
            logging.info("Load table is None, so beginning process_data")
            self.process_data()
            logging.info("Process_data done")
            logging.info("BEGIN Build table")
            self.build_table()
            logging.info("END Build table")
            logging.info("BEGIN dump data")
            self.dump_data(dump_file)
            logging.info("END dump table")
        else:
            self.load_table(load_table)

    def process_data(self):
        """
        dict_e1_e2: maps (e1,e2) to all r in data [TRAINING]
        unique_e1_r: stores unique e1_r for building table [TRAINING + TEST + VAL]
        dict_e: stores all possible e2 for a given e1 [TRAINING]
        dict_e1_r: stores all possible e2 for given e1 and r [TRAINING]
        cache_r_r1_r2: a cache for similarity(r,hadamard(r1,r2))
        """
        self.dict_e1_e2 = {}
        self.unique_e1_r = {}
        self.dict_e = {}
        #self.dict_e1_r = {}
        self.cache_r_r1_r2 = {}

        for facts in self.kb.facts:
            key = (facts[0], facts[2])
            # dict_e1_e2
            if(key not in self.dict_e1_e2):
                self.dict_e1_e2[key] = []
            self.dict_e1_e2[key].append(facts[1])

            # dict_e1_r
            # if((facts[0], facts[1]) not in self.dict_e1_r):
            #    self.dict_e1_r[(facts[0], facts[1])] = []
            #self.dict_e1_r[(facts[0], facts[1])].append(facts[2])

            # dict_e
            if facts[0] not in self.dict_e:
                self.dict_e[facts[0]] = set()
            self.dict_e[facts[0]].add(facts[2])

            # unique_e1_r being created
            if((facts[0], facts[1]) not in self.unique_e1_r):
                self.unique_e1_r[(facts[0], facts[1])] = len(self.unique_e1_r)

        for facts in self.kb_val.facts:
            if((facts[0], facts[1]) not in self.unique_e1_r):
                self.unique_e1_r[(facts[0], facts[1])] = len(self.unique_e1_r)

        for facts in self.kb_test.facts:
            if((facts[0], facts[1]) not in self.unique_e1_r):
                self.unique_e1_r[(facts[0], facts[1])] = len(self.unique_e1_r)

        #
        self.sorted_e1_r_list = list(self.unique_e1_r.keys())
        self.sorted_e1_r_list.sort()

        ####################################
        # can delete was added for debugging
        # self.r_count={}
        # for fact in self.kb.facts:
        # 	if fact[1] not in self.r_count:
        # 		self.r_count[fact[1]] = 0
        # 	self.r_count[fact[1]]+=1
        # import pdb
        # pdb.set_trace()
        ####################################

    def process(self, key):
        e1, r = key
        #global ctr
        #global total_size
        #global start_time
        # if ctr % 250 == 1:
        #    logging.info("Processed %d in %f seconds. Avg size %f" %
        #                 (ctr, time.time()-start_time, total_size/ctr))
        logging.info("Ctr: {}".format(self.ctr))
        score_dict = {}
        for u in range(len(self.kb.entity_map)):
            score, r1, u1, r2 = self.compute_score((e1, r, u))
            if(score != 0):
                score_dict[u] = (score, r1, u1, r2)

        #total_size += len(score_dict)

        if(len(score_dict.keys()) > 0):
            #self.table[(e1, r)] = score_dict
            val_list = [x[0] for x in score_dict.values()]
            mean = np.mean(val_list)
            std = np.std(val_list)
            max_score = max(val_list)
            index_max = val_list.index(max_score)
            simi_index = list(score_dict.keys())[index_max]
            stat = {"mean": mean, "std": std, "max_score": max_score,
                    "index_max": index_max, "simi_index": simi_index}
            #self.stat_table[(e1, r)] = stat
        self.ctr += 1
        return (key, score_dict, stat)

    def build_table(self):
        """
        self.table = {}
        self.stat_table = {}
        self.ctr = 0
        total_size = 0
        start_time = time.time()

        pool = multiprocessing.Pool()
        #Pdb().set_trace()
        list_of_tables = list(tqdm(pool.imap_unordered(self.process, list(self.unique_e1_r.keys())[:100]), total=100))
        #with tqdm(total= 100) as pbar:
        #    for i, _ in tqdm(enumerate(pool.imap_unordered(self.process, list(self.unique_e1_r.keys())[:100]))):
        #        pbar.update()
        #list_of_tables = pool.map(self.process,list(self.unique_e1_r.keys())[:100])
        #list_of_tables = map(self.process,list(self.unique_e1_r.keys())[:100])
        #Pdb().set_trace()
        for key, score_dict, stat in list_of_tables:
            self.table[key] = score_dict
            self.stat_table[key] = stat
        #Pdb().set_trace()
        pool.close()

        """
        # a table for each unique (e1,r)
        # Pdb().set_trace()
        entities = len(self.kb.entity_map)
        self.table = {}
        self.stat_table = {}
        ctr = 0
        start_time = time.time()
        total_size = 0

        start_index = self.offset*int(len(self.sorted_e1_r_list)/self.parts)
        end_index = (self.offset+1)*int(len(self.sorted_e1_r_list)/self.parts)
        if self.offset == self.parts-1:
            end_index = len(self.sorted_e1_r_list)

        total_e1_r = end_index - start_index
        last_logged_at = 0

        if os.path.exists(self.dump_file+'.partial'):
            logging.info("Found Partial dump at: {}. Loading from it: ".format(
                self.dump_file+'.partial'))
            self.load_table(self.dump_file+'.partial')

        logging.info("Total (e1,r): {}. Start Index: {} End Index: {}".format(
            len(self.unique_e1_r), start_index, end_index))
        for (e1, r) in tqdm(self.sorted_e1_r_list[start_index: end_index]):
            if (e1, r) in self.table:
                continue
            #
            if ctr % 250 == 1:
                logging.info("Processed %d in %f seconds. Avg size %f" %
                             (ctr, time.time()-start_time, total_size/ctr))
                start_time = time.time()
            score_dict = {}
            for u in range(entities):
                score, r1, u1, r2 = self.compute_score((e1, r, u))
                if(score != 0):
                    score_dict[u] = (score, r1, u1, r2)

            total_size += len(score_dict)

            if(len(score_dict.keys()) > 0):
                self.table[(e1, r)] = score_dict
                val_list = [x[0] for x in self.table[(e1, r)].values()]
                mean = np.mean(val_list)
                std = np.std(val_list)
                max_score = max(val_list)
                index_max = val_list.index(max_score)
                simi_index = list(self.table[(e1, r)].keys())[index_max]
                stat = {"mean": mean, "std": std, "max_score": max_score,
                        "index_max": index_max, "simi_index": simi_index}
                self.stat_table[(e1, r)] = stat
            ctr += 1
            if int(100*(ctr/total_e1_r)) >= (last_logged_at + 20):
                last_logged_at = int(100*(ctr/total_e1_r))
                logging.info(
                    "***START Creating a partial dump at : {} %".format(last_logged_at))
                self.dump_data(self.dump_file+'.partial')
                logging.info(
                    "***END Created a partial dump at : {} %".format(last_logged_at))

    def dump_data(self, filename):
        dump_dict = {}
        dump_dict['dict_e'] = self.dict_e
        dump_dict['dict_e1_e2'] = self.dict_e1_e2
        dump_dict['unique_e1_r'] = self.unique_e1_r
        #dump_dict['dict_e1_r'] = self.dict_e1_r

        dump_dict['table'] = self.table
        dump_dict['stat_table'] = self.stat_table

        with open(filename, 'wb') as inputfile:
            pickle.dump(dump_dict, inputfile)

    def load_table(self, filename):
        with open(filename, "rb") as f:
            dump_dict = pickle.load(f)
        self.dict_e = dump_dict['dict_e']
        self.dict_e1_e2 = dump_dict['dict_e1_e2']
        self.unique_e1_r = dump_dict['unique_e1_r']
        #self.dict_e1_r = dump_dict['dict_e1_r']

        self.table = dump_dict['table']
        self.stat_table = dump_dict['stat_table']

    def compute_score(self, triple):
        '''
        Returns template score for given triple
        Iterates over all e1,r depending on flag of use_hard_triple_scoring
        '''

        assert (len(triple) == 3), "Triple must contain three elements"

        best_score = 0
        best_r1 = -1
        best_u1 = -1
        best_r2 = -1

        e2 = triple[2]
        r = triple[1]
        e1 = triple[0]

        if(self.use_hard_triple_scoring == False):
            raise NotImplementedError
        else:
            entities = self.dict_e.get(e1, [])
            for u1 in entities:
                if u1 != e2:
                    relations1 = self.dict_e1_e2.get((e1, u1), [])
                    relations2 = self.dict_e1_e2.get((u1, e2), [])
                    for r1 in relations1:
                        for r2 in relations2:
                            if((r, r1, r2) in self.cache_r_r1_r2):
                                score = self.cache_r_r1_r2[(r, r1, r2)]
                            else:
                                hadamard_r1_r2 = self.base_model.get_hadamard_product(
                                    r1, r2)
                                r_emb = self.base_model.get_relation_embedding(
                                    r)
                                score = self.base_model.get_relation_similarity_from_embedding(
                                    hadamard_r1_r2, r_emb)
                                self.cache_r_r1_r2[(r, r1, r2)] = score
                            if(score > best_score):
                                best_score = score
                                best_r1 = r1
                                best_r2 = r2
                                best_u1 = u1
        return (best_score, best_r1, best_u1, best_r2)

    def get_input(self, fact):
        key = (fact[0], fact[1])
        features = [0, 0, 0, 0, 0, 0, 0]

        if(key in self.table.keys()):
            val_list = [x[0] for x in self.table[key].values()]
            if (len(val_list) != 0):
                max_score = self.stat_table[key]['max_score']
                my_score = self.table[key].get(fact[2], (0, -1, -1, -1))[0]
                simi = self.base_model.get_entity_similarity(
                    fact[2], self.stat_table[key]['simi_index'])
                rank = utils.get_rank(val_list, my_score)
                conditional_rank = rank*1.0/len(val_list)
                mean = self.stat_table[key]['mean']
                std = self.stat_table[key]['std']
                features = [my_score, max_score, simi,
                            rank, conditional_rank, mean, std]

        return features

    def get_explanation(self, fact):
        """
        returns the best two length path for this fact, score,
        best answer for (e1,r,?), best_score, best path for best answer, z_score
        we have : (score, r1, u1, r2)

        """
        key = (fact[0], fact[1])
        #features = [score, 2lp, top score, top u, top path, z score]
        features = [0,(-1,-1,-1),0,-1,(-1,-1,-1),0]
        if(key in self.table.keys()):
            my_table_stats = self.stat_table[key]
            if (len(self.table[key]) != 0):
                my_score_and_path = self.table[key].get(fact[2], (0,-1,-1,-1))
                my_score = my_score_and_path[0]
                my_best_path = my_score_and_path[1:]
                best_answer = my_table_stats['simi_index']
                best_score_and_path = self.table[key].get(best_answer)
                best_score = best_score_and_path[0]
                best_answer_path = best_score_and_path[1:]
                z_score = (my_score - my_table_stats['mean'])/(my_table_stats['std'] + utils.EPSILON)
                features = [my_score, my_best_path, best_score, best_answer,best_answer_path, z_score]

        return features

    def get_english_explanation(self, fact, explainer):
        try:
            path = self.get_explanation(fact)[1]

            explainer.check_relation(path[0],fact, source='t6r1')
            explainer.check_relation(path[2],fact, source='t6r2')
            #self.exp_template = '$html_fact_e1_r1_u1 and $html_fact_u1_r2_e2'
            return string.Template(self.exp_template).substitute(html_fact_e1_r1_u1=explainer.html_fact([fact[0], path[0], path[1]]), html_fact_u1_r2_e2=explainer.html_fact([path[1], path[2], fact[2]]))
        except:
            logging.error("Some error occured in template 6 while getting explanation\n{}".format(
                sys.exc_info()[0]))
            return explainer.NO_EXPLANATION

max_stars_repo_name: NinaCalvi/OKBC

max_issues_repo_name: NinaCalvi/OKBC

max_forks_repo_name: NinaCalvi/OKBC


text:
from __future__ import division, print_function
import numpy as np

N_GRID = 100

def isotropic_circular(R,L,nu=0.27,E=165):
    '''
    Computes the lateral stress and strain tensors for isotropic circular
    SBCA.

    Input:
    R = bending radius
    L = analyser diameter
    nu = Poisson's ratio
    E = Young's modulus

    The units of R and L do not matter, as long as they are the same. Same for
    E but the choice fixes the units of the computed stress tensor

    Output:
    X, Y = 2D matrices giving the coordinates of the surface points
    stress = The components of the stress tensor
    strain = The components of the strain tensor

    Note: stress and strain matrices contain nan for the points outside the
    crystal surface.
    '''
    #generate the grid shifts
    x=np.linspace(-L/2,L/2,N_GRID)
    X,Y=np.meshgrid(x,x)

    stress = {}
    strain = {}

    stress['xx'] = -E/(16*R**2)*(X**2 + 3*Y**2 -L**2/4)
    stress['yy'] = -E/(16*R**2)*(3*X**2 + Y**2 -L**2/4)
    stress['xy'] = E/(8*R**2)*X*Y
    stress['yx'] = stress['xy']

    stress['rr'] = E/(16*R**2)*(L**2/4-X**2-Y**2)
    stress['phiphi'] = E/(16*R**2)*(L**2/4-3*X**2-3*Y**2)
    stress['rphi'] = np.zeros(X.shape)
    stress['phir'] = stress['rphi']

    strain['xx'] = ((1-nu)*L**2/4-(1-3*nu)*X**2-(3-nu)*Y**2)/(16*R**2)
    strain['yy'] = ((1-nu)*L**2/4-(1-3*nu)*Y**2-(3-nu)*X**2)/(16*R**2)
    strain['xy'] = (1+nu)/(8*R**2)*X*Y
    strain['yx'] = strain['xy']

    strain['rr'] = 1/(16*R**2)*((1-nu)*L**2/4-(1-3*nu)*(X**2+Y**2))
    strain['phiphi'] = 1/(16*R**2)*((1-nu)*L**2/4-(3-nu)*(X**2+Y**2))
    strain['rphi'] = np.zeros(X.shape)
    strain['phir'] = strain['rphi']

    strain['zz'] = nu/(4*R**2)*(X**2+Y**2-L**2/8)

    for k in stress:
        stress[k][X**2+Y**2 > L**2/4] = np.nan
    for k in strain:
        strain[k][X**2+Y**2 > L**2/4] = np.nan

    return X,Y,stress,strain

def anisotropic_circular(R,L,S):
    '''
    Computes the lateral stress and strain tensors for anisotropic circular
    SBCA.

    Input:
    R = bending radius
    L = analyser diameter
    S = compliance matrix

    The units of R and L do not matter, as long as they are the same. Same for
    S but the choice fixes the units of the computed stress tensor

    Output:
    X, Y = 2D matrices giving the coordinates of the surface points
    stress = The components of the stress tensor
    strain = The components of the strain tensor

    Note: stress and strain matrices contain nan for the points outside the
    crystal surface.
    '''

    x=np.linspace(-L/2,L/2,N_GRID)
    X,Y=np.meshgrid(x,x)

    r_squared = X**2+Y**2
    phi = np.arctan2(Y,X)

    stress = {}
    strain = {}

    D = 1/(2*R**2*(3*(S[0,0]+S[1,1])+2*S[0,1]+S[5,5]))

    stress['xx'] = D*(L**2/4-X**2-3*Y**2)
    stress['yy'] = D*(L**2/4-3*X**2-Y**2)
    stress['xy'] = 2*D*X*Y
    stress['yx'] = stress['xy']

    stress['rr'] = D*(L**2/4-r_squared)
    stress['phiphi'] = D*(L**2/4-3*r_squared)
    stress['rphi'] = np.zeros(X.shape)
    stress['phir'] = stress['rphi']

    #shorthand notation
    uzzaux1 = (S[2,0]+S[2,1])*L**2/4
    uzzaux2 = 2*(S[2,0]+S[2,1])
    uzzaux3 = np.sqrt((S[2,1]-S[2,0])**2+S[2,5]**2)

    beta = np.arctan2(S[2,5],(S[2,1]-S[2,0]))

    strain['zz'] = D*(uzzaux1 - (uzzaux2+uzzaux3*np.sin(2*phi+beta))*r_squared)

    eff_poisson = -4*(S[2,0]+S[2,1])/(3*(S[0,0]+S[1,1])+2*S[0,1]+S[5,5])
    print('Effective Poisson ratio in lateral strain: ', str(eff_poisson))

    for k in stress:
        stress[k][X**2+Y**2 > L**2/4] = np.nan
    for k in strain:
        strain[k][X**2+Y**2 > L**2/4] = np.nan

    return X,Y,stress,strain

def isotropic_rectangular(R,a,b,nu=0.27,E=165):
    '''
    Computes the lateral stress and strain tensors for isotropic rectangular
    SBCA.

    Input:
    R = bending radius
    a, b = side lengths of the rectangle
    nu = Poisson's ratio
    E = Young's modulus

    The units of R and L do not matter, as long as they are the same. Same for
    E but the choice fixes the units of the computed stress tensor

    Output:
    X, Y = 2D matrices giving the coordinates of the surface points
           X-axis is parallel with the side a and y-axis with b
    stress = The components of the stress tensor
    strain = The components of the strain tensor

    Note: stress and strain matrices contain nan for the points outside the
    crystal surface.
    '''
    #resize grid points according to the linear dimensions so the separation
    #of the grid points is equal in x and y directions
    if a > b:
        x=np.linspace(-a/2,a/2,N_GRID)
        y=np.linspace(-b/2,b/2,int(np.round(N_GRID*b/a)))
    else:
        y=np.linspace(-b/2,b/2,N_GRID)
        x=np.linspace(-a/2,a/2,int(np.round(N_GRID*a/b)))

    X,Y=np.meshgrid(x,y)

    stress = {}
    strain = {}

    g = 8 + 10*((a/b)**2+(b/a)**2) + (1-nu)*((a/b)**2-(b/a)**2)**2

    stress['xx'] = E/(g*R**2) * (a**2/12-X**2 + ((1+nu)/2 + 5*(a/b)**2 + (1-nu)/2*(a/b)**4)*(b**2/12-Y**2))
    stress['yy'] = E/(g*R**2) * (b**2/12-Y**2 + ((1+nu)/2 + 5*(b/a)**2 + (1-nu)/2*(b/a)**2)*(a**2/12-X**2))
    stress['xy'] = 2*E/(g*R**2)*X*Y

    strain['zz'] = nu/(g*R**2) * (((3+nu)/2+5*(b/a)**2+(1-nu)/2*(b/a)**4)*(X**2 - a**2/12)+\
                                  ((3+nu)/2+5*(a/b)**2+(1-nu)/2*(a/b)**4)*(Y**2 - b**2/12))

    return X,Y,stress,strain

max_stars_repo_name: aripekka/sbcalc

max_issues_repo_name: aripekka/sbcalc

max_forks_repo_name: aripekka/sbcalc


text:
# -*- coding: utf-8 -*-
# author: Qinghe Li
# date: 2019/10/11

import os
import torch
import random
import numpy as np
from transformers.tokenization_bert import BertTokenizer

class DataLoader(object):
    def __init__(self, data_dir, bert_model_dir, params, token_pad_idx=0):
        self.data_dir = data_dir
        self.batch_size = params.batch_size
        self.max_len = params.max_len
        self.device = params.device
        self.seed = params.seed
        self.token_pad_idx = 0

        labels = self.load_labels()
        self.label2idx = {label: idx for idx, label in enumerate(labels)}
        self.idx2label = {idx: label for idx, label in enumerate(labels)}
        params.label2idx = self.label2idx
        params.idx2label = self.idx2label
        self.label_pad_idx = self.label2idx['[PAD]']

        self.tokenizer = BertTokenizer.from_pretrained(bert_model_dir, do_lower_case=True)

    def load_labels(self):
        labels = []
        file_path = os.path.join(self.data_dir, 'labels.txt')
        with open(file_path, 'r') as file:
            for label in file:
                labels.append(label.strip())
        return labels

    def load_sentences_labels(self, sentences_file, labels_file, dict):
        sentences = []
        labels = []
        with open(sentences_file, 'r', encoding='gb18030') as fs:
            lines = fs.readlines()
            for line in lines:
                tokens = line.strip().split()
                tokens.insert(0, '[CLS]')
                tokens.append('[SEP]')
                tokens_to_id =  self.tokenizer.convert_tokens_to_ids(tokens)
                sentences.append(tokens_to_id)

        with open(labels_file, 'r', encoding='gb18030') as fl:
            lines = fl.readlines()
            for line in lines:
                label = line.strip().split()
                label.insert(0, '[CLS]')
                label.append('[SEP]')
                label_to_id = [self.label2idx[key] for key in label]
                labels.append(label_to_id)

        # checks to ensure there is a label for each token
        assert len(sentences) == len(labels)
        for i in range(len(sentences)):
            assert len(labels[i]) == len(sentences[i])

        # storing sentences and labels in dict d
        dict['data'] = sentences
        dict['labels'] = labels
        dict['size'] = len(sentences)

    def load_data(self, data_type):
        data = {}

        sentences_file = os.path.join(self.data_dir, data_type, 'sentences.txt')
        labels_path = os.path.join(self.data_dir, data_type, 'labels.txt')
        self.load_sentences_labels(sentences_file, labels_path, data)

        return data

    def data_iterator(self, data, shuffle=False):
        # make a list that decides the order in which we go over the data- this avoids explicit shuffling of data
        order = list(range(data['size']))
        if shuffle:
            random.seed(self.seed)
            random.shuffle(order)

        # one pass over data
        for i in range(data['size'] // self.batch_size):
            # fetch sentences and labels
            sentences = [data['data'][idx] for idx in order[i*self.batch_size: (i+1)*self.batch_size]]
            labels = [data['labels'][idx] for idx in order[i*self.batch_size: (i+1)*self.batch_size]]

            # batch length
            batch_len = len(sentences)

            # compute length of longest sentence in batch
            batch_max_len = max([len(s) for s in sentences])
            max_len = min(batch_max_len, self.max_len)

            # prepare a numpy array with the data, initialising the data with pad_idx
            batch_data = self.token_pad_idx * np.ones((batch_len, max_len))
            batch_labels = self.label_pad_idx * np.ones((batch_len, max_len))

            # copy the data to the numpy array
            for j in range(batch_len):
                cur_len = len(sentences[j])
                if cur_len <= max_len:
                    batch_data[j][: cur_len] = sentences[j]
                    batch_labels[j][: cur_len] = labels[j]
                else:
                    batch_data[j] = sentences[j][: max_len]
                    batch_labels[j] = labels[j][: max_len]

            # since all data are indices, we convert them to torch LongTensors
            batch_data = torch.tensor(batch_data, dtype=torch.long)
            batch_labels = torch.tensor(batch_labels, dtype=torch.long)

            # shift tensors to GPU if available
            batch_data, batch_labels = batch_data.to(self.device), batch_labels.to(self.device)

            yield batch_data, batch_labels

max_stars_repo_name: Derrors/Others

max_issues_repo_name: Derrors/Others

max_forks_repo_name: Derrors/Others