text:
r"""

Autoregressive (AR) processes are processes of the form:

.. math::

  x(n) = a(1)x(n-1) + a(2)x(n-2) + ... + a(P)x(n-P) + e(n)

where e(n) is a white noise process. The usage of 'e' suggests interpreting
the linear combination of P past values of x(n) as the minimum mean square
error linear predictor of x(n) Thus

.. math::

  e(n) = x(n) - a(1)x(n-1) - a(2)x(n-2) - ... - a(P)x(n-P)

Due to whiteness, e(n) is also pointwise uncorrelated--ie,

.. math::
   :nowrap:

   \begin{align*}
      \text{(i)}   && E\{e(n)e^{*}(n-m)\}& = \delta(n-m) &\\
      \text{(ii)}  && E\{e(n)x^{*}(m)\}  & = 0           & m\neq n\\
      \text{(iii)} && E\{|e|^{2}\} = E\{e(n)e^{*}(n)\} &= E\{e(n)x^{*}(n)\} &
   \end{align*}

These principles form the basis of the methods in this module for
estimating the AR coefficients and the error/innovations power.
"""

import numpy as np
from nitime.lazy import scipy_linalg as linalg

import nitime.utils as utils
from spectral import freq_response

def AR_est_YW(x, order, rxx=None):
    r"""Determine the autoregressive (AR) model of a random process x using
    the Yule Walker equations. The AR model takes this convention:

    .. math::

      x(n) = a(1)x(n-1) + a(2)x(n-2) + \dots + a(p)x(n-p) + e(n)

    where e(n) is a zero-mean white noise process with variance sig_sq,
    and p is the order of the AR model. This method returns the a_i and
    sigma

    The orthogonality property of minimum mean square error estimates
    states that

    .. math::

      E\{e(n)x^{*}(n-k)\} = 0 \quad 1\leq k\leq p

    Inserting the definition of the error signal into the equations above
    yields the Yule Walker system of equations:

    .. math::

      R_{xx}(k) = \sum_{i=1}^{p}a(i)R_{xx}(k-i) \quad1\leq k\leq p

    Similarly, the variance of the error process is

    .. math::

      E\{e(n)e^{*}(n)\}   = E\{e(n)x^{*}(n)\} = R_{xx}(0)-\sum_{i=1}^{p}a(i)R^{*}(i)

    Parameters
    ----------
    x : ndarray
        The sampled autoregressive random process

    order : int
        The order p of the AR system

    rxx : ndarray (optional)
        An optional, possibly unbiased estimate of the autocorrelation of x

    Returns
    -------
    ak, sig_sq: The estimated AR coefficients and innovations variance

    """
    if rxx is not None and type(rxx) == np.ndarray:
        r_m = rxx[:order + 1]
    else:
        r_m = utils.autocorr(x)[:order + 1]

    Tm = linalg.toeplitz(r_m[:order])
    y = r_m[1:]
    ak = linalg.solve(Tm, y)
    sigma_v = r_m[0].real - np.dot(r_m[1:].conj(), ak).real
    return ak, sigma_v

def AR_est_LD(x, order, rxx=None):
    r"""Levinson-Durbin algorithm for solving the Hermitian Toeplitz
    system of Yule-Walker equations in the AR estimation problem

    .. math::

       T^{(p)}a^{(p)} = \gamma^{(p+1)}

    where

    .. math::
       :nowrap:

       \begin{align*}
       T^{(p)} &= \begin{pmatrix}
          R_{0} & R_{1}^{*} & \cdots & R_{p-1}^{*}\\
          R_{1} & R_{0} & \cdots & R_{p-2}^{*}\\
          \vdots & \vdots & \ddots & \vdots\\
          R_{p-1}^{*} & R_{p-2}^{*} & \cdots & R_{0}
       \end{pmatrix}\\
       a^{(p)} &=\begin{pmatrix} a_1 & a_2 & \cdots a_p \end{pmatrix}^{T}\\
       \gamma^{(p+1)}&=\begin{pmatrix}R_1 & R_2 & \cdots & R_p \end{pmatrix}^{T}
       \end{align*}

    and :math:`R_k` is the autocorrelation of the kth lag

    Parameters
    ----------

    x: ndarray
      the zero-mean stochastic process
    order : int
      the AR model order--IE the rank of the system.
    rxx : ndarray, optional
      (at least) order+1 samples of the autocorrelation sequence

    Returns
    -------

    ak, sig_sq
      The AR coefficients for 1 <= k <= p, and the variance of the
      driving white noise process

    """

    if rxx is not None and type(rxx) == np.ndarray:
        rxx_m = rxx[:order + 1]
    else:
        rxx_m = utils.autocorr(x)[:order + 1]
    w = np.zeros((order + 1, ), rxx_m.dtype)
    # intialize the recursion with the R[0]w[1]=r[1] solution (p=1)
    b = rxx_m[0].real
    w_k = rxx_m[1] / b
    w[1] = w_k
    p = 2
    while p <= order:
        b *= 1 - (w_k * w_k.conj()).real
        w_k = (rxx_m[p] - (w[1:p] * rxx_m[1:p][::-1]).sum()) / b
        # update w_k from k=1,2,...,p-1
        # with a correction from w*_i i=p-1,p-2,...,1
        w[1:p] = w[1:p] - w_k * w[1:p][::-1].conj()
        w[p] = w_k
        p += 1
    b *= 1 - (w_k * w_k.conj()).real
    return w[1:], b

def lwr_recursion(r):
    r"""Perform a Levinson-Wiggins[Whittle]-Robinson recursion to
    find the coefficients a(i) that satisfy the matrix version
    of the Yule-Walker system of P + 1 equations:

    sum_{i=0}^{P} a(i)r(k-i) = 0, for k = {1,2,...,P}

    with the additional equation

    sum_{i=0}^{P} a(i)r(-k) = V

    where V is the covariance matrix of the innovations process,
    and a(0) is fixed at the identity matrix

    Also note that r is defined as:

    r(k) = E{ X(t)X*(t-k) } ( * = conjugate transpose )
    r(-k) = r*(k)

    This routine adapts the algorithm found in eqs (1)-(11)
    in Morf, Vieira, Kailath 1978

    Parameters
    ----------

    r : ndarray, shape (P + 1, nc, nc)

    Returns
    -------

    a : ndarray (P,nc,nc)
      coefficient sequence of order P
    sigma : ndarray (nc,nc)
      covariance estimate

    """

    # r is (P+1, nc, nc)
    nc = r.shape[1]
    P = r.shape[0] - 1

    a = np.zeros((P, nc, nc))  # ar coefs
    b = np.zeros_like(a)  # lp coefs
    sigb = np.zeros_like(r[0])  # forward prediction error covariance
    sigf = np.zeros_like(r[0])  # backward prediction error covariance
    delta = np.zeros_like(r[0])

    # initialize
    idnt = np.eye(nc)
    sigf[:] = r[0]
    sigb[:] = r[0]

    # iteratively find sequences A_{p+1}(i) and B_{p+1}(i)
    for p in xrange(P):

        # calculate delta_{p+1}
        # delta_{p+1} = r(p+1) + sum_{i=1}^{p} a(i)r(p+1-i)
        delta[:] = r[p + 1]
        for i in xrange(1, p + 1):
            delta += np.dot(a[i - 1], r[p + 1 - i])

        # intermediate values XXX: should turn these into solution-problems
        ka = np.dot(delta, linalg.inv(sigb))
        kb = np.dot(delta.conj().T, linalg.inv(sigf))

        # store a_{p} before updating sequence to a_{p+1}
        ao = a.copy()
        # a_{p+1}(i) = a_{p}(i) - ka*b_{p}(p+1-i) for i in {1,2,...,p}
        # b_{p+1}(i) = b_{p}(i) - kb*a_{p}(p+1-i) for i in {1,2,...,p}
        for i in xrange(1, p + 1):
            a[i - 1] -= np.dot(ka, b[p - i])
        for i in xrange(1, p + 1):
            b[i - 1] -= np.dot(kb, ao[p - i])

        a[p] = -ka
        b[p] = -kb

        sigf = np.dot(idnt - np.dot(ka, kb), sigf)
        sigb = np.dot(idnt - np.dot(kb, ka), sigb)

    return a, sigf

def MAR_est_LWR(x, order, rxx=None):
    r"""
    MAR estimation, using the LWR algorithm, as in Morf et al.

    Parameters
    ----------
    x : ndarray
        The sampled autoregressive random process

    order : int
        The order P of the AR system

    rxx : ndarray (optional)
        An optional, possibly unbiased estimate of the autocovariance of x

    Returns
    -------
    a, ecov: The system coefficients and the estimated covariance
    """
    Rxx = utils.autocov_vector(x, nlags=order)
    a, ecov = lwr_recursion(Rxx.transpose(2, 0, 1))
    return a, ecov

def AR_psd(ak, sigma_v, n_freqs=1024, sides='onesided'):
    r"""
    Compute the PSD of an AR process, based on the process coefficients and
    covariance

    n_freqs : int
        The number of spacings on the frequency grid from [-PI,PI).
        If sides=='onesided', n_freqs/2+1 frequencies are computed from [0,PI]

    sides : str (optional)
        Indicates whether to return a one-sided or two-sided PSD

    Returns
    -------
    (w, ar_psd)
    w : Array of normalized frequences from [-.5, .5) or [0,.5]
    ar_psd : A PSD estimate computed by sigma_v / |1-a(f)|**2 , where
             a(f) = DTFT(ak)

    """
    # compute the psd as |H(f)|**2, where H(f) is the transfer function
    # for this model s[n] = a1*s[n-1] + a2*s[n-2] + ... aP*s[n-P] + v[n]
    # Taken as a IIR system with unit-variance white noise input e[n]
    # and output s[n],
    # b0*e[n] = w0*s[n] + w1*s[n-1] + w2*s[n-2] + ... + wP*s[n-P],
    # where b0 = sqrt(VAR{v[n]}), w0 = 1, and wk = -ak for k>0
    # the transfer function here is H(f) = DTFT(w)
    # leading to Sxx(f)/Exx(f) = |H(f)|**2 = VAR{v[n]} / |W(f)|**2
    w, hw = freq_response(sigma_v ** 0.5, a=np.r_[1, -ak],
                          n_freqs=n_freqs, sides=sides)
    ar_psd = (hw * hw.conj()).real
    return (w, 2 * ar_psd) if sides == 'onesided' else (w, ar_psd)

#-----------------------------------------------------------------------------
# Granger causality analysis
#-----------------------------------------------------------------------------
def transfer_function_xy(a, n_freqs=1024):
    r"""Helper routine to compute the transfer function H(w) based
    on sequence of coefficient matrices A(i). The z transforms
    follow from this definition:

    X[t] + sum_{k=1}^P a[k]X[t-k] = Err[t]

    Parameters
    ----------

    a : ndarray, shape (P, 2, 2)
      sequence of coef matrices describing an mAR process
    n_freqs : int, optional
      number of frequencies to compute in range [0,PI]

    Returns
    -------

    Hw : ndarray
      The transfer function from innovations process vector to
      mAR process X

    """
    # these concatenations follow from the observation that A(0) is
    # implicitly the identity matrix
    ai = np.r_[1, a[:, 0, 0]]
    bi = np.r_[0, a[:, 0, 1]]
    ci = np.r_[0, a[:, 1, 0]]
    di = np.r_[1, a[:, 1, 1]]

    # compute A(w) such that A(w)X(w) = Err(w)
    w, aw = freq_response(ai, n_freqs=n_freqs)
    _, bw = freq_response(bi, n_freqs=n_freqs)
    _, cw = freq_response(ci, n_freqs=n_freqs)
    _, dw = freq_response(di, n_freqs=n_freqs)

    #A = np.array([ [1-aw, -bw], [-cw, 1-dw] ])
    A = np.array([[aw, bw], [cw, dw]])
    # compute the transfer function from Err to X. Since Err(w) is 1(w),
    # the transfer function H(w) = A^(-1)(w)
    # (use 2x2 matrix shortcut)
    detA = (A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0])
    Hw = np.array([[dw, -bw], [-cw, aw]])
    Hw /= detA
    return w, Hw

def spectral_matrix_xy(Hw, cov):
    r"""Compute the spectral matrix S(w), from the convention:

    X[t] + sum_{k=1}^P a[k]X[t-k] = Err[t]

    The formulation follows from Ding, Chen, Bressler 2008,
    pg 6 eqs (11) to (15)

    The transfer function H(w) should be computed first from
    transfer_function_xy()

    Parameters
    ----------

    Hw : ndarray (2, 2, n_freqs)
      Pre-computed transfer function from transfer_function_xy()

    cov : ndarray (2, 2)
      The covariance between innovations processes in Err[t]

    Returns
    -------

    Sw: ndarrays
      matrix of spectral density functions
    """

    nw = Hw.shape[-1]
    # now compute specral density function estimate
    # S(w) = H(w)SigH*(w)
    Sw = np.empty((2, 2, nw), 'D')

    # do a shortcut for 2x2:
    # compute T(w) = SigH*(w)
    # t00 = Sig[0,0] * H*_00(w) + Sig[0,1] * H*_10(w)
    t00 = cov[0, 0] * Hw[0, 0].conj() + cov[0, 1] * Hw[0, 1].conj()
    # t01 = Sig[0,0] * H*_01(w) + Sig[0,1] * H*_11(w)
    t01 = cov[0, 0] * Hw[1, 0].conj() + cov[0, 1] * Hw[1, 1].conj()
    # t10 = Sig[1,0] * H*_00(w) + Sig[1,1] * H*_10(w)
    t10 = cov[1, 0] * Hw[0, 0].conj() + cov[1, 1] * Hw[0, 1].conj()
    # t11 = Sig[1,0] * H*_01(w) + Sig[1,1] * H*_11(w)
    t11 = cov[1, 0] * Hw[1, 0].conj() + cov[1, 1] * Hw[1, 1].conj()

    # now S(w) = H(w)T(w)
    Sw[0, 0] = Hw[0, 0] * t00 + Hw[0, 1] * t10
    Sw[0, 1] = Hw[0, 0] * t01 + Hw[0, 1] * t11
    Sw[1, 0] = Hw[1, 0] * t00 + Hw[1, 1] * t10
    Sw[1, 1] = Hw[1, 0] * t01 + Hw[1, 1] * t11

    return Sw

def coherence_from_spectral(Sw):
    r"""Compute the spectral coherence between processes X and Y,
    given their spectral matrix S(w)

    Parameters
    ----------

    Sw : ndarray
      spectral matrix
    """

    Sxx = Sw[0, 0].real
    Syy = Sw[1, 1].real

    Sxy_mod_sq = (Sw[0, 1] * Sw[1, 0]).real
    Sxy_mod_sq /= Sxx
    Sxy_mod_sq /= Syy
    return Sxy_mod_sq

def interdependence_xy(Sw):
    r"""Compute the 'total interdependence' between processes X and Y,
    given their spectral matrix S(w)

    Parameters
    ----------

    Sw : ndarray
      spectral matrix

    Returns
    -------

    fxy(w)
      interdependence function of frequency
    """

    Cw = coherence_from_spectral(Sw)
    return -np.log(1 - Cw)

def granger_causality_xy(a, cov, n_freqs=1024):
    r"""Compute the Granger causality between processes X and Y, which
    are linked in a multivariate autoregressive (mAR) model parameterized
    by coefficient matrices a(i) and the innovations covariance matrix

    X[t] + sum_{k=1}^P a[k]X[t-k] = Err[t]

    Parameters
    ----------

    a : ndarray, (P,2,2)
      coefficient matrices characterizing the autoregressive mixing
    cov : ndarray, (2,2)
      covariance matrix characterizing the innovations vector
    n_freqs: int
      number of frequencies to compute in the fourier transform

    Returns
    -------

    w, f_x_on_y, f_y_on_x, f_xy, Sw
      1) vector of frequencies
      2) function of the Granger causality of X on Y
      3) function of the Granger causality of Y on X
      4) function of the 'instantaneous causality' between X and Y
      5) spectral density matrix
    """

    w, Hw = transfer_function_xy(a, n_freqs=n_freqs)

    sigma = cov[0, 0]
    upsilon = cov[0, 1]
    gamma = cov[1, 1]

    # this transformation of the transfer functions computes the
    # Granger causality of Y on X
    gamma2 = gamma - upsilon ** 2 / sigma

    Hxy = Hw[0, 1]
    Hxx_hat = Hw[0, 0] + (upsilon / sigma) * Hxy

    xx_auto_component = (sigma * Hxx_hat * Hxx_hat.conj()).real
    cross_component = gamma2 * Hxy * Hxy.conj()
    Sxx = xx_auto_component + cross_component
    f_y_on_x = np.log(Sxx.real / xx_auto_component)

    # this transformation computes the Granger causality of X on Y
    sigma2 = sigma - upsilon ** 2 / gamma

    Hyx = Hw[1, 0]
    Hyy_hat = Hw[1, 1] + (upsilon / gamma) * Hyx
    yy_auto_component = (gamma * Hyy_hat * Hyy_hat.conj()).real
    cross_component = sigma2 * Hyx * Hyx.conj()
    Syy = yy_auto_component + cross_component
    f_x_on_y = np.log(Syy.real / yy_auto_component)

    # now compute cross densities, using the latest transformation
    Hxx = Hw[0, 0]
    Hyx = Hw[1, 0]
    Hxy_hat = Hw[0, 1] + (upsilon / gamma) * Hxx
    Sxy = sigma2 * Hxx * Hyx.conj() + gamma * Hxy_hat * Hyy_hat.conj()
    Syx = sigma2 * Hyx * Hxx.conj() + gamma * Hyy_hat * Hxy_hat.conj()

    # can safely throw away imaginary part
    # since Sxx and Syy are real, and Sxy == Syx*
    detS = (Sxx * Syy - Sxy * Syx).real
    f_xy = xx_auto_component * yy_auto_component
    f_xy /= detS
    f_xy = np.log(f_xy)

    return w, f_x_on_y, f_y_on_x, f_xy, np.array([[Sxx, Sxy], [Syx, Syy]])

max_stars_repo_name: mluessi/nitime

max_issues_repo_name: mluessi/nitime

max_forks_repo_name: mluessi/nitime


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

# Import modules

import sys
import os
import glob
import xarray as xa
import xesmf as xe
import pandas as pd
from itertools import product
# import dask.array as da
import numpy as np
import re
from dask.distributed import Client
from rcat.utils import ini_reader
from rcat.utils.polygons import mask_region
import rcat.runtime.RCAT_stats as st
import rcat.utils.grids as gr

import warnings
warnings.filterwarnings("ignore")

# Functions
def get_settings(config_file):
    """
    Retrieve information from main configuration file
    """
    conf_dict = ini_reader.get_config_dict(config_file)
    d = {
        'models': conf_dict['MODELS'],
        'obs metadata file': conf_dict['OBS']['metadata file'],
        'obs start year': conf_dict['OBS']['start year'],
        'obs end year': conf_dict['OBS']['end year'],
        'obs months': conf_dict['OBS']['months'],
        'variables': conf_dict['SETTINGS']['variables'],
        'var modification': conf_dict['SETTINGS']['variable modification'],
        'regions': conf_dict['SETTINGS']['regions'],
        'requested_stats': conf_dict['STATISTICS']['stats'],
        'stats_conf': st.mod_stats_config(conf_dict['STATISTICS']['stats']),
        'validation plot': conf_dict['PLOTTING']['validation plot'],
        'map configure': conf_dict['PLOTTING']['map configure'],
        'map grid setup': conf_dict['PLOTTING']['map grid setup'],
        'map kwargs': conf_dict['PLOTTING']['map kwargs'],
        'line grid setup': conf_dict['PLOTTING']['line grid setup'],
        'line kwargs': conf_dict['PLOTTING']['line kwargs'],
        'cluster type': conf_dict['CLUSTER']['cluster type'],
        'nodes': conf_dict['CLUSTER']['nodes'],
        'cluster kwargs': conf_dict['CLUSTER']['cluster kwargs'],
        'outdir': conf_dict['SETTINGS']['output dir'],
    }

    return d

def local_cluster_setup():
    """
    Set up local-pc cluster

    """
    from dask.distributed import LocalCluster
    cluster = LocalCluster(processes=False)
    return cluster

def slurm_cluster_setup(nodes=1, **kwargs):
    """
    Set up SLURM cluster

    Parameters
    ----------
    nodes: int
        Number of nodes to use
    **kwargs:
        Keyword arguments for cluster specifications
    """
    from dask_jobqueue import SLURMCluster
    cluster = SLURMCluster(**kwargs)
    cluster.scale(nodes)
    return cluster

def get_args():
    """
    Read configuration file
    Parameters
    ----------
    -
    Returns
    -------
    Input arguments
    """
    import argparse

    # Configuring argument setup and handling
    parser = argparse.ArgumentParser(
        description='Main script for model/obs validation')
    parser.add_argument('--config', '-c',  metavar='name config file',
                        type=str, help='<Required> Full path to config file',
                        required=True)
    return parser.parse_args()

def get_grid_coords(nc, grid_coords):
    """
    Read model grid coordinates

    Parameters
    ----------
    nc: xarray dataset
    Returns
    -------
    grid_coords: dict
    """

    def _domain(lats, lons, x):
        lons_p = np.r_[lons[x, x:-x], lons[x:-x, -1-x][1:],
                       lons[-1-x, x:-x][::-1][1:], lons[x:-x, x][::-1][1:]]
        lats_p = np.r_[lats[x, x:-x], lats[x:-x, -1-x][1:],
                       lats[-1-x, x:-x][::-1][1:], lats[x:-x, x][::-1][1:]]
        return list(zip(lons_p, lats_p))

    # If lon/lat is 1D; create 2D meshgrid
    lons = nc.lon.values
    lats = nc.lat.values
    lon, lat = np.meshgrid(lons, lats)\
        if lats.ndim == 1 else (lons, lats)

    # Calculate domain mid point if not given
    idx = tuple([np.int(i/2) for i in lat.shape])
    lat0 = lat[idx]
    lon0 = lon[idx]

    grid_coords['lat_0'] = lat0
    grid_coords['lon_0'] = lon0

    gp_bfr = 1
    grid_coords['crnrs'] = [lat[gp_bfr, gp_bfr], lon[gp_bfr, gp_bfr],
                            lat[-gp_bfr, -gp_bfr], lon[-gp_bfr, -gp_bfr]]
    grid_coords['domain'] = _domain(lat, lon, gp_bfr)

    return grid_coords

def get_grids(nc, target_grid, method='bilinear'):
    """
    Get and/or modify the source and target grids for interpolation
    """
    # EDIT 15/11 2019 by Petter Lind
    # There's no consistent way to handle data where lon_bnds/lat_bnds are
    # already available, They may have different formats not all handled by
    # xesmf regridding tool. Thus, for now, calculation of grid corners is
    # always performed when method='conservative'
    if method == 'conservative':
        slon_b, slat_b = gr.fnCellCorners(nc.lon.values, nc.lat.values)
        s_grid = {'lon': nc.lon.values, 'lat': nc.lat.values,
                  'lon_b': slon_b, 'lat_b': slat_b}
        tlon_b, tlat_b = gr.fnCellCorners(target_grid['lon'],
                                          target_grid['lat'])
        t_grid = {'lon': target_grid['lon'], 'lat': target_grid['lat'],
                  'lon_b': tlon_b, 'lat_b': tlat_b}
    else:
        s_grid = {'lon': nc.lon.values, 'lat': nc.lat.values}
        t_grid = target_grid

    return s_grid, t_grid

def regrid_func(dd, v, vconf, mnames, onames, gdict):
    if vconf['regrid'] is None:
        gdict.update(
            {'lon': {m: dd[m]['grid']['lon'] for m in mnames},
             'lat': {m: dd[m]['grid']['lat'] for m in mnames}})
        if None not in onames:
            for obs in onames:
                gdict['lon'].update({obs: dd[obs]['grid']['lon']})
                gdict['lat'].update({obs: dd[obs]['grid']['lat']})
        gridname = 'native_grid'
    else:
        if isinstance(vconf['regrid'], dict):
            target_grid = xa.open_dataset(vconf['regrid']['file'])
            gridname = "grid_{}".format(vconf['regrid']['name'])
            gdict.update({'lon': {gridname: target_grid['lon'].values},
                          'lat': {gridname: target_grid['lat'].values}})
            for mod in mnames:
                dd[mod]['data'] = regrid_calc(dd, mod, v, target_grid,
                                              vconf['rgr method'])
            if None not in onames:
                for obs in onames:
                    dd[obs]['data'] = regrid_calc(dd, obs, v, target_grid,
                                                  vconf['rgr method'])
        elif vconf['regrid'] in onames:
            oname = vconf['regrid']
            target_grid = dd[oname]['grid']
            gridname = dd[oname]['gridname']
            gdict.update({'lon': {oname: target_grid['lon']},
                          'lat': {oname: target_grid['lat']}})
            for m in mnames:
                dd[m]['data'] = regrid_calc(dd, m, v, target_grid,
                                            vconf['rgr method'])
            if len(onames) > 1:
                obslist = onames.copy()
                obslist.remove(oname)
                for obs in obslist:
                    dd[obs]['data'] = regrid_calc(dd, obs, v, target_grid,
                                                  vconf['rgr method'])
        elif vconf['regrid'] in mnames:
            mname = vconf['regrid']
            modlist = mnames.copy()
            modlist.remove(mname)
            target_grid = dd[mname]['grid']
            gridname = dd[mname]['gridname']
            gdict.update({'lon': {mname: target_grid['lon']},
                          'lat': {mname: target_grid['lat']}})
            for mod in modlist:
                dd[mod]['data'] = regrid_calc(dd, mod, v, target_grid,
                                              vconf['rgr method'])
            if None not in onames:
                for obs in onames:
                    dd[obs]['data'] = regrid_calc(dd, obs, v, target_grid,
                                                  vconf['rgr method'])
        else:
            raise ValueError(("\n\n\tTarget grid name not found!\n"
                              "Check 'regrid to' option in main config file"))

    gdict.update({'gridname': gridname})

    return dd, gdict

def regrid_calc(data, data_name, var, target_grid, method):
    print("\t\t** Regridding {} data **\n".format(data_name.upper()))

    indata = data[data_name]['data']

    # Get grid info
    sgrid, tgrid = get_grids(indata, target_grid, method)

    # Regridding
    regridder = gr.add_matrix_NaNs(xe.Regridder(sgrid, tgrid, method))
    # regridder.clean_weight_file()
    ds_rgr = regridder(indata[var])
    dset = ds_rgr.to_dataset()

    return dset

def resampling(data, v, tresample):
    """
    Resample data to chosen time frequency and resample method.
    """
    # if isinstance(tresample, dict):
    #     for tr, val in tresample.items():
    #         lval = val if isinstance(val, list) else [val]
    #         if tr == 'select hours':
    #             data = data.sel(time=np.isin(data['time.hour'], lval))
    #         elif tr == 'select dates':
    #             data = data.sel(time=lval)
    #         else:
    #             errmsg = (f"\t\t\nUnknown type of selection: {tr}!\n Only "
    #                      f"'select hours' and 'select dates' are available.")
    #             raise ValueError(errmsg)
    errmsg = (f"\t\t\nResample resolution argument must be a list of two "
              f"items. Please re-check the configuration.")
    assert len(tresample) == 2, errmsg

    if tresample[0] in ('select hours', 'select dates'):
        val = tresample[1]
        lval = val if isinstance(val, list) else [val]
        if tresample[0] == 'select hours':
            data = data.sel(time=np.isin(data['time.hour'], lval))
        else:
            data = data.sel(time=lval)
    else:
        from pandas import to_timedelta
        diff = data.time.values[1] - data.time.values[0]
        nsec = to_timedelta(diff).total_seconds()
        tr, fr = _get_freq(tresample[0])
        sec_resample = to_timedelta(tr, fr).total_seconds()
        if nsec != sec_resample:
            data = eval("data.resample(time='{}').{}('time').\
                          dropna('time', 'all')".format(
                              tresample[0], tresample[1]))
        else:
            print("\t\tData is already at target resolution, skipping "
                  "resampling ...\n")
    return data

def conditional_data(dd_condition, cond_var, cond_data, indata):
    """
    Sub-select data conditional on other data
    """
    import operator

    ops = {'>': operator.gt,
           '<': operator.lt,
           '>=': operator.ge,
           '<=': operator.le,
           '=': operator.eq}

    def _get_cond_mask(darr, thr_type, relate, q, expand=None, axis=0):
        if thr_type == 'percentile':
            mask = np.apply_along_axis(
                lambda x: ops[relate](x, np.percentile(x, q=q)),
                axis=axis, arr=darr)
        elif thr_type == 'static':
            mask = np.apply_along_axis(lambda x: ops[relate](x, q),
                                       axis=axis, arr=darr)
        if expand is not None:
            mask = np.repeat(mask, expand, axis=axis)

        return mask

    # Resample data
    if 'resample resolution' in dd_condition:
        print("\t\tResampling conditional data ...\n")
        tres = dd_condition['resample resolution']
        cond_data = resampling(cond_data, cond_var, tres)

    # Re-chunk in space dim if needed
    if len(cond_data.chunks['time']) > 1:
        cond_data = manage_chunks(cond_data, 'space')

    # Type of conditional
    cond_type = dd_condition['type']

    # Logical operator for conditional selection
    relate = dd_condition['operator']

    if cond_type == 'file':
        errmsg = ("\t\t\nConditional data from file must be 2D!\n")
        with xa.open_dataset(dd_condition['file in']) as fopen:
            q = fopen[dd_condition['file var']].squeeze()
            assert q.ndim == 2, errmsg
        sub_data = indata.where(ops[relate](indata, q))
    else:
        q = float(dd_condition['value'])

        if cond_data.time.size != indata.time.size:
            expand = indata.time.size/cond_data.time.size
        else:
            expand = None

        _mask = xa.apply_ufunc(
            _get_cond_mask, cond_data[cond_var], input_core_dims=[['time']],
            output_core_dims=[['time']], dask='parallelized',
            output_dtypes=[float], exclude_dims={'time'},
            dask_gufunc_kwargs={'output_sizes': {'time': indata.time.size}},
            kwargs={'thr_type': cond_type, 'relate': relate, 'q': q,
                    'expand': expand, 'axis': -1})
        dims = list(_mask.dims)
        mask = _mask.transpose(dims[-1], dims[0], dims[1])
        sub_data = indata.where(mask)

    return sub_data

def calc_stats(ddict, vlist, stat, pool, chunk_dim, stats_config, regions):
    """
    Calculate statistics for variables and models/obs
    """
    st_data = {}
    for v in vlist:
        st_data[v] = {}
        for m in ddict[v]:
            print("\tCalculate {} {} for {}\n".format(v, stat, m))

            if ddict[v][m] is None:
                print("\t* Data not available for model {}, moving on".format(
                    m))
            else:
                indata = ddict[v][m]['data']
                st_data[v][m] = {}
                if len(indata.data_vars) > 2:
                    indata = indata[v].to_dataset()

                # Resampling of data
                if stats_config[stat]['resample resolution'] is not None:
                    resample_args = stats_config[stat]['resample resolution']
                    if isinstance(resample_args, dict):
                        if v in resample_args:
                            print("\t\tResampling input data ...\n")
                            indata = resampling(indata, v, resample_args[v])
                        else:
                            pass
                    else:
                        print("\t\tResampling input data ...\n")
                        indata = resampling(indata, v, resample_args)

                # Check chunking of data
                data = manage_chunks(indata, chunk_dim)
                # print(f"Chunks of {m}: {data.chunks}")

                # Conditional analysis; sample data according to condition
                # (static threshold,  percentile or from file).
                if stats_config[stat]['cond analysis'] is not None:
                    cond_dict = stats_config[stat]['cond analysis']
                    if v in cond_dict:
                        print("\t\tPerforming conditional sub-selection\n")
                        cond_calc = cond_dict[v]
                        cond_var = cond_calc['cvar']

                        if cond_var == v:
                            cond_data = indata
                        else:
                            # assert vcond in st_data, msg.format(vcond, v)
                            cond_data = ddict[cond_var][m]['data']

                        # Extract sub selection of data
                        sub_data = conditional_data(cond_calc, cond_var,
                                                    cond_data, indata)

                        data = manage_chunks(sub_data, chunk_dim)

                # Calculate stats
                st_data[v][m]['domain'] = st.calc_statistics(data, v, stat,
                                                             stats_config)

                if regions:
                    masks = {
                        r: mask_region(
                            data.lon.values, data.lat.values, r,
                            cut_data=False) for r in regions}
                    if pool:
                        mdata = {r: get_masked_data(data, v, masks[r])
                                 for r in regions}

                        st_mdata = {r: st.calc_statistics(mdata[r], v, stat,
                                                          stats_config)
                                    for r in regions}
                    else:
                        st_mdata = {r: get_masked_data(st_data[v][m]['domain'],
                                                       v, masks[r])
                                    for r in regions}

                    st_data[v][m]['regions'] = st_mdata

    return st_data

def get_month_string(mlist):
    d = {1: 'J', 2: 'F', 3: 'M', 4: 'A', 5: 'M', 6: 'J', 7: 'J', 8: 'A',
         9: 'S', 10: 'O', 11: 'N', 12: 'D'}
    dl = {1: 'Jan', 2: 'Feb', 3: 'Mar', 4: 'Apr', 5: 'May', 6: 'Jun', 7: 'Jul',
          8: 'Aug', 9: 'Sep', 10: 'Oct', 11: 'Nov', 12: 'Dec'}
    if len(mlist) == 12:
        mstring = 'ANN'
    elif len(mlist) == 1:
        mstring = dl[mlist[0]]
    else:
        if mlist in ([1, 2, 12], (1, 2, 12)):
            mlist = (12, 1, 2)
        mstring = ''.join(d[m] for m in mlist)
    return mstring

def _get_freq(tf):
    from functools import reduce

    d = [j.isdigit() for j in tf]
    if np.any(d):
        freq = int(reduce((lambda x, y: x+y), [x for x, y in zip(tf, d) if y]))
    else:
        freq = 1

    unit = reduce((lambda x, y: x+y), [x for x, y in zip(tf, d) if not y])

    if unit in ('M', 'Y'):
        freq = freq*30 if unit == 'M' else freq*365
        unit = 'D'
    elif unit[0] == 'Q':
        # Quarterly frequencies - assuming average 90 days
        freq = 90
        unit = 'D'

    return freq, unit

def _space_dim(ds):
    """
    Return labels for space dimensions in data set.

    Space dimensions in different observations have different names. This
    dictionary is created to account for that (to some extent), but in the
    future this might be change/removed. For now, it's hard coded.
    """
    spcdim = {'x': ['x', 'X', 'lon', 'rlon'],
              'y': ['y', 'Y', 'lat', 'rlat']}
    xd = [x for x in ds.dims if x in spcdim['x']][0]
    yd = [y for y in ds.dims if y in spcdim['y']][0]

    return xd, yd

def save_to_disk(data, label, stat, odir, var, grid, sy, ey, tsuffix,
                 stat_dict, tres, thr='', regs=None):
    """
    Saving statistical data to netcdf files
    """
    # encoding = {'lat': {'dtype': 'float32', '_FillValue': False},
    #             'lon': {'dtype': 'float32', '_FillValue': False},
    #             var: {'dtype': 'float32', '_FillValue': 1.e20}
    #             }
    if stat in ('annual cycle', 'seasonal cycle', 'diurnal cycle'):
        # _tstat = stat_dict['stat method'].partition(' ')[0]
        # tstat = '_' + re.sub(r'[^a-zA-Z0-9 \.\n]', '', _tstat).\
        #         replace(' ', '_')
        tstat = '_' + stat_dict['stat method'].replace(' ', '')
    else:
        tstat = ''
    stat_name = stat.replace(' ', '_')
    if stat in ('diurnal cycle'):
        stat_fn = "{}_{}".format(stat_name, stat_dict['dcycle stat'])
    else:
        stat_fn = stat_name

    if regs is not None:
        for r in regs:
            rn = r.replace(' ', '_')
            data['regions'][r].attrs['Analysed time'] =\
                "{}-{} | {}".format(sy, ey, tsuffix)
            fname = '{}_{}_{}_{}{}{}_{}_{}_{}-{}_{}.nc'.format(
                label, stat_fn, var, thr, tres, tstat, rn, grid, sy, ey,
                tsuffix)
            data['regions'][r].to_netcdf(os.path.join(odir, stat_name, fname))
    fname = '{}_{}_{}_{}{}{}_{}_{}-{}_{}.nc'.format(
        label, stat_fn, var, thr, tres, tstat, grid, sy, ey, tsuffix)
    data['domain'].attrs['Analysed time'] = "{}-{} | {}".format(sy, ey,
                                                                tsuffix)
    data['domain'].to_netcdf(os.path.join(odir, stat_name, fname))

def get_masked_data(data, var, mask):
    """
    Mask region
    """
    mask_in = xa.DataArray(np.broadcast_to(mask, data[var].shape),
                           dims=data[var].dims)
    return data.where(mask_in, drop=True)

    # Petter Lind 2019-10-29
    # N.B. For masking of large data sets, the code below might be a more
    # viable choice.
    #
    # def _mask_func(arr, axis=0, lons=None, lats=None, region=''):
    #     iter_3d = arr.shape[axis]
    #     mdata = mask_region(lons, lats, region, arr, iter_3d=iter_3d,
    #                               cut_data=True)[0]
    #     return mdata

    # mask = mask_region(data.lon.values, data.lat.values, reg,
    #                          cut_data=True)

    # imask = mask[0]
    # mask_data = da.map_blocks(_mask_func, data[var].data, dtype=float,
    #                           chunks=(data.chunks['time'],
    #                                   imask.shape[0], imask.shape[1]),
    #                           lons=data.lon.values, lats=data.lat.values,
    #                           region=reg)
    # lon_d = 'x' if mask[1].ndim == 1 else ['y', 'x']
    # lat_d = 'y' if mask[2].ndim == 1 else ['y', 'x']
    # out = xa.Dataset(
    #     {var: (['time', 'y', 'x'],  mask_data)},
    #     coords={'lon': (lon_d, mask[1]), 'lat': (lat_d, mask[2]),
    #             'time': data.time.values})
    # return out

def manage_chunks(data, chunk_dim):
    """
    Re-chunk the data in specified dimension.
    """
    xd, yd = _space_dim(data)
    xsize = data[xd].size
    ysize = data[yd].size

    # Rule of thumb: chunk size should at least be 1e6 elements
    # http://xarray.pydata.org/en/stable/dask.html#chunking-and-performance
    # Max chunksize just an arbitrary value to assert 'reasonable' chunking
    min_chunksize = 1e6
    max_chunksize = 1e8
    if chunk_dim == 'space':
        chunksize = np.mean(data.chunks[xd]) * np.mean(data.chunks[yd])\
                * data.time.size
        if chunksize < min_chunksize or chunksize > max_chunksize:
            sub_size = np.sqrt(min_chunksize/data.time.size)
            _csize_x = int((xsize/sub_size))
            csize_x = 5 if _csize_x < 5 else _csize_x
            _csize_y = int((ysize/sub_size))
            csize_y = 5 if _csize_y < 5 else _csize_y
            data = data.chunk({'time': -1, xd: csize_x, yd: csize_y})
        else:
            data = data.chunk({'time': -1})
    else:
        chunksize = xsize * ysize * np.mean(data.chunks['time'])
        if chunksize < min_chunksize or chunksize > max_chunksize:
            sub_size = max_chunksize/(xsize*ysize)
            _csize_t = int(data.time.size/sub_size)
            csize_t = 5 if _csize_t < 5 else _csize_t
            data = data.chunk({'time': csize_t, xd: xsize, yd: ysize})
        else:
            data = data.chunk({xd: xsize, yd: ysize})

    return data

def get_variable_config(var_config, var):
    """
    Retrieve configuration info for variable var as defined in main
    configuration file,
    """
    vdict = {
        'var names': var_config['var names'],
        'input resolution': var_config['freq'],
        'units': var_config['units'],
        'scale factor': var_config['scale factor'],
        'obs scale factor': var_config['obs scale factor'] if
        'obs scale factor' in var_config else None,
        'deacc': var_config['accumulated'],
        'regrid': var_config['regrid to'] if 'regrid to' in
        var_config else None,
        'rgr method': var_config['regrid method'] if 'regrid method' in
        var_config else None,
    }
    return vdict

def variabel_modification(dd, nv_dd, funargs, mlist, olist):
    """
    Create new or modify existing variables.
    """
    from inspect import signature
    out_dd = {}
    expression = nv_dd['expression']
    func = eval(f"lambda {funargs}: {expression}")
    args = list(signature(func).parameters)
    _modlist = nv_dd['models']
    _obslist = nv_dd['obs']
    if _modlist is not None:
        modlist = _modlist if isinstance(_modlist, list) else mlist \
                if _modlist == 'all' else [_modlist]
        for m in modlist:
            input_data = [dd[nv_dd['input'][a]][m]['data'][nv_dd['input'][a]]
                          for a in args]
            _new_data = xa.apply_ufunc(
                func, *input_data, input_core_dims=[[]]*len(input_data),
                dask='parallelized', output_dtypes=[float],
            )
            new_data = _new_data.to_dataset(name=new_var)
            out_dd[m] = {'data': new_data}
    if _obslist is not None:
        obslist = _obslist if isinstance(_obslist, list) else olist \
                if _obslist == 'all' else [_obslist]
        for o in obslist:
            input_data = {k: dd[v][o]['data'][v]
                          for k, v in nv_dd['input'].items()}
            _new_data = func(**input_data)
            new_data = _new_data.to_dataset(name=new_var)
            out_dd[o] = {'data': new_data}

    return out_dd

def get_mod_data(model, mconf, tres, var, vnames, cfactor, deacc):
    """
    Open model data files where file path is dependent on time resolution tres.
    """
    import re

    print("\t-- Opening {} files\n".format(model.upper()))

    fyear = mconf['start year']
    lyear = mconf['end year']
    months = mconf['months']
    date_list = ["{}{:02d}".format(yy, mm) for yy, mm in product(
        range(fyear, lyear+1), months)]

    if vnames is not None:
        if 'all' in vnames:
            readvar = vnames['all']['prfx']
        elif model in vnames:
            readvar = vnames[model]['prfx']
        else:
            readvar = var
    else:
        readvar = var
    file_path = os.path.join(mconf['fpath'],
                             f'{tres}/{readvar}/{readvar}_*.nc')
    _flist = glob.glob(file_path)

    errmsg = (f"Could not find any files at specified location:\n{file_path} "
              "\n\nexiting ...")
    if not _flist:
        print("\t\n{}".format(errmsg))
        sys.exit()

    _file_dates = [re.split('-|_', f.rsplit('.')[-2])[-2:] for f in _flist]
    file_dates = [(d[0][:6], d[1][:6]) for d in _file_dates]
    fidx = [np.where([d[0] <= date <= d[1] for d in file_dates])[0]
            for date in date_list]
    flist = [_flist[i] for i in np.unique(np.hstack(fidx))]
    flist.sort()

    if np.unique([len(f) for f in flist]).size > 1:
        flngth = np.unique([len(f) for f in flist])
        flist = [f for f in flist if len(f) == flngth[0]]

    # Chunk sizes in temporal and spatial dimensions
    ch_t = mconf['chunks_time']
    ch_x = mconf['chunks_x']
    ch_y = mconf['chunks_y']

    if deacc:
        _mdata = xa.open_mfdataset(
            flist, parallel=True,  # engine='h5netcdf',
            data_vars='minimal', coords='minimal',
            concat_dim='time', combine='by_coords', compat='override',
            chunks={**ch_t, **ch_x, **ch_y},
            preprocess=(lambda arr: arr.diff('time')))
    else:
        _mdata = xa.open_mfdataset(
            flist, parallel=True,  # engine='h5netcdf',
            data_vars='minimal', coords='minimal',
            concat_dim='time', combine='by_coords', compat='override',
            chunks={**ch_t, **ch_x, **ch_y})

    # Time stamps
    if 'units' in _mdata.time.attrs:
        if _mdata.time.attrs['units'] == 'day as %Y%m%d.%f':
            dates = [pd.to_datetime(d, format='%Y%m%d') +
                     pd.to_timedelta((d % 1)*86400, unit='s').round('H')
                     for d in _mdata.time.values]
            _mdata['time'] = (('time',), dates)
        else:
            raise ValueError(("\n\n\tUnfortunately, at the moment the time "
                              "units in file cannot be treated, change if "
                              "possible"))

    # Extract years
    mdata = _mdata.where(((_mdata.time.dt.year >= fyear) &
                          (_mdata.time.dt.year <= lyear) &
                          (np.isin(_mdata.time.dt.month, months))), drop=True)

    # Remove height dim
    if 'height' in mdata.dims:
        mdata = mdata.squeeze()

    # Rename variable if not consistent with name in configuration file
    if vnames is not None:
        if 'all' in vnames:
            mdata = mdata.rename({vnames['all']['vname']: var})
        elif model in vnames:
            mdata = mdata.rename({vnames[model]['vname']: var})

    if cfactor is not None:
        mdata[var] *= cfactor

    # Model grid
    gridname = 'grid_{}'.format(mconf['grid name'])

    # - Unrotate grid if needed
    if mconf['grid type'] == 'rot':
        lon_reg, lat_reg = gr.rotated_grid_transform(
            mdata.rlon, mdata.rlat,
            mdata.rotated_pole.grid_north_pole_longitude,
            mdata.rotated_pole.grid_north_pole_latitude)
        mdata = mdata.assign_coords({'lon': (('y', 'x'), lon_reg),
                                     'lat': (('y', 'x'), lat_reg)}).\
            swap_dims({'rlon': 'x', 'rlat': 'y'})
        grid = {'lon': lon_reg, 'lat': lat_reg}
    else:
        grid = {'lon': mdata.lon.values, 'lat': mdata.lat.values}

    outdata = {'data': mdata.unify_chunks(),
               'grid': grid, 'gridname': gridname}

    return outdata

def get_obs_data(metadata_file, obs, var, cfactor, sy, ey, mns):
    """
    Open obs data files.
    """
    from importlib.machinery import SourceFileLoader
    obs_meta = SourceFileLoader("obs_meta", metadata_file).load_module()
    # obs_dict = obs_meta.obs_data()

    sdate = '{}{:02d}'.format(sy, np.min(mns))
    edate = '{}{:02d}'.format(ey, np.max(mns))

    print("\t-- Opening {} files\n".format(obs.upper()))

    # Open obs files
    obs_flist = obs_meta.get_file_list(var, obs, sdate, edate)

    emsg = ("Could not find any {} files at specified location"
            "\n\nexiting ...".format(obs.upper()))
    if not obs_flist:
        print("\t\n{}".format(emsg))
        sys.exit()

    date_list = ["{}{:02d}".format(yy, mm) for yy, mm in product(
        range(sy, ey+1), mns)]
    _file_dates = [re.split('-|_', f.rsplit('.')[-2])[-2:] for f in obs_flist]
    file_dates = [(d[0][:6], d[1][:6]) for d in _file_dates]
    fidx = [np.where([d[0] <= date <= d[1] for d in file_dates])[0]
            for date in date_list]
    flist = [obs_flist[i] for i in np.unique(np.hstack(fidx))]
    flist.sort()
    f_obs = xa.open_mfdataset(
        flist, parallel=True,   # engine='h5netcdf',
        data_vars='minimal', coords='minimal',
        concat_dim='time', combine='by_coords',
        compat='override').unify_chunks()

    # Extract years
    obs_data = f_obs.where(((f_obs.time.dt.year >= sy) &
                            (f_obs.time.dt.year <= ey) &
                            (np.isin(f_obs.time.dt.month, mns))),
                           drop=True)
    # Scale data
    if cfactor is not None:
        obs_data[var] *= cfactor

    lons = obs_data.lon.values
    lats = obs_data.lat.values

    # Drop bnds dims
    if 'bnds' in obs_data.dims:
        obs_data = obs_data.drop_dims('bnds')

    # Labels for space dimensions
    xd, yd = _space_dim(obs_data)

    # Make sure lon/lat elements are in ascending order
    if lons.ndim == 1:
        if np.diff(lons)[0] < 0:
            lons = lons[::-1]
            obs_data = obs_data.reindex({xd: obs_data[xd][::-1]})
        if np.diff(lats)[0] < 0:
            lats = lats[::-1]
            obs_data = obs_data.reindex({yd: obs_data[yd][::-1]})
    elif lons.ndim == 2:
        if np.diff(lons[0, :])[0] < 0:
            lons = np.flipud(lons)
            obs_data = obs_data.reindex({xd: np.flipud(obs_data[xd])})
        if np.diff(lats[:, 0])[0] < 0:
            lats = np.flipud(lats)
            obs_data = obs_data.reindex({yd: np.flipud(obs_data[yd])})

    grid = {'lon': lons, 'lat': lats}
    gridname = 'grid_{}'.format(obs.upper())

    outdata = {'data': obs_data, 'grid': grid, 'gridname': gridname}
    return outdata

def get_plot_dict(cdict, var, grid_coords, models, obs, yrs_d, mon_d, tres,
                  img_outdir, stat_outdir, stat):
    """
    Create a dictionary with settings for a validation plot procedure.
    """
    # Get some settings and meta data
    vconf = get_variable_config(cdict['variables'][var], var)
    grdnme = grid_coords['target grid'][var]['gridname']
    st = stat.replace(' ', '_')
    if stat == 'diurnal cycle':
        stnm = "{}_{}".format(st, cdict['stats_conf'][stat]['dcycle stat'])
    else:
        stnm = st
    if stat in ('annual cycle', 'seasonal cycle', 'diurnal cycle'):
        # if 'percentile' in cdict['stats_conf'][stat]['stat method']:
        #     tstat = '_percentile'
        tstat = '_' + cdict['stats_conf'][stat]['stat method'].replace(
            ' ', '')
    else:
        tstat = ''
    thrlg = (('thr' in cdict['stats_conf'][stat]) and
             (cdict['stats_conf'][stat]['thr'] is not None) and
             (var in cdict['stats_conf'][stat]['thr']))
    if thrlg:
        thrstr = "thr{}_".format(str(cdict['stats_conf'][stat]['thr'][v]))
    else:
        thrstr = ''

    # Create dictionaries with list of files for models and obs
    _fm_list = {stat: [glob.glob(os.path.join(
        stat_outdir, '{}'.format(st), '{}_{}_{}_{}{}{}_{}_{}-{}_{}.nc'.format(
            m, stnm, var, thrstr, tres, tstat, grdnme, yrs_d[m][0],
            yrs_d[m][1], get_month_string(mon_d[m])))) for m in models]}
    fm_list = {s: [y for x in ll for y in x] for s, ll in _fm_list.items()}

    obs_list = [obs] if not isinstance(obs, list) else obs
    if obs is not None:
        _fo_list = {stat: [glob.glob(os.path.join(
            stat_outdir, '{}'.format(st),
            '{}_{}_{}_{}{}{}_{}_{}-{}_{}.nc'.format(
                o, stnm, var, thrstr, tres, tstat, grdnme,
                yrs_d[o][0], yrs_d[o][1], get_month_string(mon_d[o]))))
            for o in obs_list]}
        fo_list = {s: [y for x in ll for y in x] for s, ll in _fo_list.items()}
    else:
        fo_list = {stat: None}

    sy = np.unique([n['start year'] for m, n in cdict['models'].items()])
    ey = np.unique([n['end year'] for m, n in cdict['models'].items()])
    if sy.size > 1:
        years = {'T1': (sy[0], ey[0]), 'T2': (sy[1], ey[1])}
    else:
        years = (sy[0], ey[0])

    # Create a dictionary with settings for plotting procedure
    plot_dict = {
        'mod files': fm_list,
        'obs files': fo_list,
        'models': models,
        'observation': obs,
        'variable': var,
        'time res': tres,
        'stat method': tstat,
        'units': vconf['units'],
        'grid coords': grid_coords,
        'map configure': cdict['map configure'],
        'map grid setup': cdict['map grid setup'],
        'map kwargs': cdict['map kwargs'],
        'line grid setup': cdict['line grid setup'],
        'line kwargs': cdict['line kwargs'],
        'regions': cdict['regions'],
        'years': years,
        'img dir': os.path.join(img_outdir, st)
    }

    # If there are regions, create list of files for these  as well
    # Then also update plot dictionary
    if cdict['regions'] is not None:
        _fm_listr = {stat: {r:  [glob.glob(os.path.join(
            stat_outdir, '{}'.format(st),
            '{}_{}_{}_{}{}{}_{}_{}_{}-{}_{}.nc'.format(
                m, stnm, var, thrstr, tres, tstat, r.replace(' ', '_'), grdnme,
                yrs_d[m][0], yrs_d[m][1], get_month_string(mon_d[m]))))
            for m in models] for r in cdict['regions']}}
        fm_listr = {s: {r: [y for x in _fm_listr[s][r] for y in x]
                        for r in _fm_listr[s]} for s in _fm_listr}
        if obs is not None:
            _fo_listr = {stat: {r: [glob.glob(os.path.join(
                stat_outdir, '{}'.format(st),
                '{}_{}_{}_{}{}{}_{}_{}_{}-{}_{}.nc'.format(
                    o, stnm, var, thrstr, tres, tstat, r.replace(' ', '_'),
                    grdnme, yrs_d[o][0], yrs_d[o][1],
                    get_month_string(mon_d[o])))) for o in obs_list]
                for r in cdict['regions']}}
            fo_listr = {s: {r: [y for x in _fo_listr[s][r] for y in x]
                            for r in _fo_listr[s]} for s in _fo_listr}
        else:
            fo_listr = {stat: None}

        plot_dict.update({
            'mod reg files': fm_listr,
            'obs reg files': fo_listr,
        })

    return plot_dict

# Configuration
args = get_args()

# Read config file
config_file = args.config
if not os.path.isfile(config_file):
    raise ValueError(f"\nConfig file, '{config_file}', does not exist!")
cdict = get_settings(config_file)

###############################################
#                                             #
#            0. CREATE DIRECTORIES            #
#                                             #
###############################################

# Create dirs
stat_outdir = os.path.join(cdict['outdir'], 'stats')
stat_names = [s.replace(' ', '_') for s in cdict['requested_stats']]
img_outdir = os.path.join(cdict['outdir'], 'imgs')

if os.path.exists(cdict['outdir']):
    msg = ("\nOutput folder\n\n{}\n\nalready exists!\nDo you want "
           "to overwrite? y/n: ".format(cdict['outdir']))
    overwrite = input(msg)
    if overwrite == 'y':
        [os.makedirs(os.path.join(stat_outdir, t), exist_ok=True)
         for t in stat_names]
        [os.makedirs(os.path.join(img_outdir, t), exist_ok=True)
         for t in stat_names]
    else:
        sys.exit()
else:
    [os.makedirs(os.path.join(stat_outdir, t)) for t in stat_names]
    [os.makedirs(os.path.join(img_outdir, t)) for t in stat_names]

# Set up distributed client
if cdict['cluster type'] == 'local':
    cluster = local_cluster_setup()
    client = Client(cluster)
elif cdict['cluster type'] == 'slurm':
    nnodes = cdict['nodes']
    sl_kwargs = cdict['cluster kwargs']
    cluster = slurm_cluster_setup(nodes=nnodes, **sl_kwargs)
    client = Client(cluster)
else:
    print("\n\tCluster type not implemented! Exiting..")
    sys.exit()

###############################################
#                                             #
#       1. OPEN DATA FILES AND REMAPPING      #
#                                             #
###############################################

# Loop over time resolution of input data
grid_coords = {}
grid_coords['meta data'] = {}
grid_coords['target grid'] = {}
data_dict = {}
month_dd = {}
year_dd = {}

print("\n=== READ & PRE-PROCESS DATA ===")
for var in cdict['variables']:
    data_dict[var] = {}
    var_conf = get_variable_config(cdict['variables'][var], var)

    tres = cdict['variables'][var]['freq']
    print("\n\tvariable: {}  |  input resolution: {}".format(var, tres))

    grid_coords['meta data'][var] = {}
    grid_coords['target grid'][var] = {}

    obs_metadata_file = cdict['obs metadata file']
    obs_name = cdict['variables'][var]['obs']
    obs_list = [obs_name] if not isinstance(obs_name, list) else obs_name
    obs_scf = var_conf['obs scale factor']
    obs_scf = ([obs_scf]*len(obs_list)
               if not isinstance(obs_scf, list) else obs_scf)
    if obs_name is not None:
        for oname, cfactor in zip(obs_list, obs_scf):
            obs_data = get_obs_data(
                obs_metadata_file, oname, var, cfactor,
                cdict['obs start year'],
                cdict['obs end year'], cdict['obs months'])
            data_dict[var][oname] = obs_data

    mod_names = []
    mod_scf = var_conf['scale factor']
    mod_scf = ([mod_scf]*len(cdict['models'].keys())
               if not isinstance(mod_scf, list) else mod_scf)
    for (mod_name, settings), scf in zip(cdict['models'].items(), mod_scf):
        mod_names.append(mod_name)
        mod_data = get_mod_data(
            mod_name, settings, tres, var, var_conf['var names'],
            scf, var_conf['deacc'])
        data_dict[var][mod_name] = mod_data

        # Update grid information for plotting purposes
        if mod_name not in grid_coords['meta data'][var]:
            grid_coords['meta data'][var][mod_name] = {}
            get_grid_coords(mod_data['data'],
                            grid_coords['meta data'][var][mod_name])

    month_dd[var] = {m: cdict['models'][m]['months'] for m in mod_names}
    year_dd[var] = {m: (cdict['models'][m]['start year'],
                        cdict['models'][m]['end year']) for m in mod_names}
    if obs_list is not None:
        month_dd[var].update({o: cdict['obs months'] for o in obs_list})
        year_dd[var].update({o: (cdict['obs start year'],
                                 cdict['obs end year']) for o in obs_list})

    ##################################
    #  1B REMAP DATA TO COMMON GRID  #
    ##################################
    regrid_func(data_dict[var], var, var_conf, mod_names, obs_list,
                grid_coords['target grid'][var])

#############################################
#  1C MODIFICATION & CREATION OF VARIABLES  #
#############################################
if cdict['var modification'] is not None:
    # Modification/manipulation of variables
    for new_var, nv_dict in cdict['var modification'].items():
        arglist = list(nv_dict['input'].keys())
        inargs = ",".join(arglist)
        data_dict[new_var] = variabel_modification(data_dict, nv_dict, inargs,
                                                   mod_names, obs_list)

        # Change parameters and dictionaries accordingly
        cdict['variables'][new_var] = \
            cdict['variables'][nv_dict['input'][arglist[0]]]
        month_dd[new_var] = month_dd[nv_dict['input'][arglist[0]]]
        year_dd[new_var] = year_dd[nv_dict['input'][arglist[0]]]
        grid_coords['target grid'][new_var] =\
            grid_coords['target grid'][nv_dict['input'][arglist[0]]]
        grid_coords['meta data'][new_var] =\
            grid_coords['meta data'][nv_dict['input'][arglist[0]]]

        # Remove input variables if 'replace' in config_main.ini is set to True
        if nv_dict['replace']:
            [data_dict.pop(v) for k, v in nv_dict['input'].items()]
            [cdict['variables'].pop(v) for k, v in nv_dict['input'].items()]

###############################################
#                                             #
#          2. CALCULATE STATISTICS            #
#                                             #
###############################################

print("\n=== STATISTICS ===\n")
stats_dict = {}
for stat in cdict['stats_conf']:
    vrs = cdict['stats_conf'][stat]['vars']
    varlist = list(cdict['variables'].keys()) if not vrs \
        else [vrs] if not isinstance(vrs, list) else vrs
    pool = cdict['stats_conf'][stat]['pool data']
    chunk_dim = cdict['stats_conf'][stat]['chunk dimension']
    stats_dict[stat] = calc_stats(data_dict, varlist, stat, pool,
                                  chunk_dim, cdict['stats_conf'],
                                  cdict['regions'])

###############################################
#                                             #
#           3. WRITE DATA TO DISK             #
#                                             #
###############################################

print("\n=== SAVE OUTPUT ===")
tresstr = {}
for stat in cdict['stats_conf']:
    tresstr[stat] = {}
    resample_res = cdict['stats_conf'][stat]['resample resolution']
    for v in stats_dict[stat]:
        thrlg = (('thr' in cdict['stats_conf'][stat]) and
                 (cdict['stats_conf'][stat]['thr'] is not None) and
                 (v in cdict['stats_conf'][stat]['thr']))
        if thrlg:
            thrval = str(cdict['stats_conf'][stat]['thr'][v])
            thrstr = "thr{}_".format(thrval)
        else:
            thrstr = ''
        if resample_res is not None:
            if isinstance(resample_res, dict):
                if v in resample_res:
                    if resample_res[v][0] in ('select dates', 'select hours'):
                        tresstr[stat][v] = resample_res[v][0].replace(' ', '_')
                    else:
                        tresstr[stat][v] = "_".join(resample_res[v])
                else:
                    tresstr[stat][v] = cdict['variables'][v]['freq']
            else:
                if resample_res[0] in ('select dates', 'select hours'):
                    tresstr[stat][v] = resample_res[0].replace(' ', '_')
                else:
                    tresstr[stat][v] = "_".join(resample_res)
        else:
            tresstr[stat][v] = cdict['variables'][v]['freq']
        gridname = grid_coords['target grid'][v]['gridname']

        for m in stats_dict[stat][v]:
            print(f"\n\twriting {m.upper()} - {v} - {stat} to disk ...")
            time_suffix = get_month_string(month_dd[v][m])
            st_data = stats_dict[stat][v][m]

            save_to_disk(st_data, m, stat, stat_outdir, v, gridname,
                         year_dd[v][m][0], year_dd[v][m][1], time_suffix,
                         cdict['stats_conf'][stat], tresstr[stat][v], thrstr,
                         cdict['regions'])

###############################################
#                                             #
#                 4. PLOTTING                 #
#                                             #
###############################################

if cdict['validation plot']:
    print('\n=== PLOTTING ===')
    import rcat.runtime.RCAT_plots as rplot
    statnames = list(stats_dict.keys())
    for sn in statnames:
        for v in stats_dict[sn]:
            plot_dict = get_plot_dict(cdict, v, grid_coords, mod_names,
                                      cdict['variables'][v]['obs'], year_dd[v],
                                      month_dd[v], tresstr[sn][v], img_outdir,
                                      stat_outdir, sn)
            print("\t\n** Plotting: {} for {} **".format(sn, v))
            rplot.plot_main(plot_dict, sn)

ext: py

max_stars_repo_name: HCLIMcom/RCAT

max_issues_repo_name: HCLIMcom/RCAT

max_forks_repo_name: HCLIMcom/RCAT

reason: import numpy


text:
import streamlit as st
import pandas as pd
import numpy as np
import joblib
from sklearn.ensemble import ExtraTreesClassifier
from prediction import get_prediction, ordinal_encoder

model = joblib.load(r'Notebook/Extratreeclassifier_model.joblib')

st.set_page_config(page_title="Accident Severity Prediction App",
                   page_icon="🚧", layout="wide")

#creating option list for dropdown menu
options_day = ['Sunday', "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday"]
options_age = ['18-30', '31-50', 'Over 51', 'Unknown', 'Under 18']

options_acc_area = ['Other', 'Office areas', 'Residential areas', ' Church areas',
       ' Industrial areas', 'School areas', '  Recreational areas',
       ' Outside rural areas', ' Hospital areas', '  Market areas',
       'Rural village areas', 'Unknown', 'Rural village areasOffice areas',
       'Recreational areas']

options_cause = ['No distancing', 'Changing lane to the right',
       'Changing lane to the left', 'Driving carelessly',
       'No priority to vehicle', 'Moving Backward',
       'No priority to pedestrian', 'Other', 'Overtaking',
       'Driving under the influence of drugs', 'Driving to the left',
       'Getting off the vehicle improperly', 'Driving at high speed',
       'Overturning', 'Turnover', 'Overspeed', 'Overloading', 'Drunk driving',
       'Unknown', 'Improper parking']
options_vehicle_type = ['Automobile', 'Lorry (41-100Q)', 'Other', 'Pick up upto 10Q',
       'Public (12 seats)', 'Stationwagen', 'Lorry (11-40Q)',
       'Public (13-45 seats)', 'Public (> 45 seats)', 'Long lorry', 'Taxi',
       'Motorcycle', 'Special vehicle', 'Ridden horse', 'Turbo', 'Bajaj', 'Bicycle']
options_driver_exp = ['5-10yr', '2-5yr', 'Above 10yr', '1-2yr', 'Below 1yr', 'No Licence', 'unknown']
options_lanes = ['Two-way (divided with broken lines road marking)', 'Undivided Two way',
       'other', 'Double carriageway (median)', 'One way',
       'Two-way (divided with solid lines road marking)', 'Unknown']

features = ['hour','day_of_week','casualties','accident_cause','vehicles_involved','vehicle_type','driver_age','accident_area','driving_experience','lanes']

st.markdown("
Accident Severity Prediction App 🚧
", unsafe_allow_html=True)
def main():
    with st.form('prediction_form'):

        st.subheader("Enter the input for following features:")

        hour = st.slider("Pickup Hour: ", 0, 23, value=0, format="%d")
        day_of_week = st.selectbox("Select Day of the Week: ", options=options_day)
        casualties = st.slider("Hour of Accident: ", 1, 8, value=0, format="%d")
        accident_cause = st.selectbox("Select Accident Cause: ", options=options_cause)
        vehicles_involved = st.slider("Pickup Hour: ", 1, 7, value=0, format="%d")
        vehicle_type = st.selectbox("Select Vehicle Type: ", options=options_vehicle_type)
        driver_age = st.selectbox("Select Driver Age: ", options=options_age)
        accident_area = st.selectbox("Select Accident Area: ", options=options_acc_area)
        driving_experience = st.selectbox("Select Driving Experience: ", options=options_driver_exp)
        lanes = st.selectbox("Select Lanes: ", options=options_lanes)

        submit = st.form_submit_button("Predict")

    if submit:
        day_of_week = ordinal_encoder(day_of_week, options_day)
        accident_cause = ordinal_encoder(accident_cause, options_cause)
        vehicle_type = ordinal_encoder(vehicle_type, options_vehicle_type)
        driver_age =  ordinal_encoder(driver_age, options_age)
        accident_area =  ordinal_encoder(accident_area, options_acc_area)
        driving_experience = ordinal_encoder(driving_experience, options_driver_exp)
        lanes = ordinal_encoder(lanes, options_lanes)

        data = np.array([hour,day_of_week,casualties,accident_cause,vehicles_involved,
                            vehicle_type,driver_age,accident_area,driving_experience,lanes]).reshape(1,-1)

        pred = get_prediction(data =data, model=model)

        st.write(f"The predicted severity is:  {pred[0]}")

if __name__ == '__main__':
    main()

max_stars_repo_name: KarthickTk/RTA-project

max_issues_repo_name: KarthickTk/RTA-project

max_forks_repo_name: KarthickTk/RTA-project


text:
import sys
import googlemaps
import pyowm
import pandas as pd
import numpy as np
import geopy.distance
from selenium import webdriver
from selenium.webdriver.chrome.options import Options
from selenium.webdriver.support.ui import Select
import os
import time
import datetime
from flask import Flask, render_template, request, jsonify, abort

OWM_API_KEY = sys.argv[1] #weather api key should be passed as first argument
MAP_API_KEY = sys.argv[2] #gmap api key should be passed as second argument

owm = pyowm.OWM(OWM_API_KEY)
gmaps = googlemaps.Client(key=MAP_API_KEY)
stops = pd.DataFrame()

# instantiate a chrome options object so you can set the size and headless preference
chrome_options = Options()
chrome_options.add_argument("--headless")
chrome_options.add_argument("--window-size=1920x1080")

app = Flask(__name__)

@app.route('/')
def serve_app():
    return render_template('roger.html')

@app.route('/roger/api/combined')
def serve_combined():
    lat = request.args.get('lat')
    lng = request.args.get('lng')
    if(lat==None or lng==None):
        abort(422)
    else:
        weather = getWeather(float(lat),float(lng))
        bus = scrape_bus_data(float(lat),float(lng))
        response = jsonify({'weather':weather,'bus':bus})
        response.headers.add('Access-Control-Allow-Origin', '*')
        return response

@app.route('/roger/api/weather')
def serve_weather_data():
    lat = request.args.get('lat')
    lng = request.args.get('lng')
    if(lat==None or lng==None):
        abort(422)
    else:
        status = getWeather(float(lat),float(lng))
        return jsonify({'status': status})

@app.route('/roger/api/bus')
def serve_bus_data():
    lat = request.args.get('lat')
    lng = request.args.get('lng')
    if(lat==None or lng==None):
        abort(422)
    else:
        bus_info = scrape_bus_data(float(lat),float(lng))
        return jsonify({'bus_info': bus_info})

'''Finds the nearest bus stop to the lat lng passed in and returns a dictioanry
containing the nearest stop name and another dictionary of atmost 5 buses arriving
at the bus stop'''
def scrape_bus_data(lat,lng):
    cls_stop_code = closest_stop(lat,lng)
    cls_stop_name = stops.loc[cls_stop_code]['stop_name']
    print(cls_stop_name)

    '''Get the buses at this stop'''
    # chromedriver binary MUST be put in the same directory as this file
    chrome_driver = os.getcwd()+'/chromedriver';

    # go to transit website
    driver = webdriver.Chrome(chrome_options=chrome_options, executable_path=chrome_driver)
    driver.get("https://truetime.portauthority.org/bustime/eta/eta.jsp")

    #find stop number input box
    stpNumInput = driver.find_element_by_id("txtFindStopID")

    #enter stop number
    stpNumInput.send_keys(cls_stop_code)

    #click the find button
    findStop_button = driver.find_element_by_id("findStop")
    findStop_button.click()

    #get the first bus data
    bus_data = {}
    for i in range(1,6):
        route = "route" + str(i)
        eta = "time" + str(i)
        busNumber = (driver.find_element_by_id(route).text).strip()
        eta = (driver.find_element_by_id(eta).text).strip()
        if busNumber != '' and eta != '':
            bus_data[i] = {'bus':busNumber, 'eta':eta}

    #return bus and stop Data
    data = {}
    data['stp_name'] = cls_stop_name
    data['bus_data'] = bus_data
    return data

'''Finds the closet stop to the given lat, lng'''
def closest_stop(lat,lng):
    global stops
    """ Loop through the stops to find the closest """
    #read stop data only once
    if stops.empty:
        getStopInfo()

    cl_stp = None #closet stop
    min = 99999999999999 # a very big value for distance
    myloc = (float(lat),float(lng))

    for index, row in stops.iterrows():
        lt = row['stop_lat']
        ln = row['stop_lon']
        stoploc = (lt,ln)
        distance = geopy.distance.vincenty(myloc, stoploc).miles
        if(distance<min):
            min = distance
            cl_stp = index

    print(cl_stp)
    return cl_stp

'''Calls the weather api to get weather at a particular latitute and
longitutde '''
def getWeather(lat,lng):
    obs = owm.weather_at_coords(lat,lng)
    #print(lat,',',lng,':',obs.get_weather().get_status())
    #return obs.get_weather().get_status()
    obs_weather = obs.get_weather()
    weather = {}
    weather['status'] = obs_weather.get_status()
    weather['detailed_status'] = obs_weather.get_detailed_status()
    weather['weather_code'] = obs_weather.get_weather_code()
    weather['weather_icon_name'] = obs_weather.get_weather_icon_name()
    weather['icon'] = 'http://openweathermap.org/img/w/'+weather['weather_icon_name']+'.png'
    weather['recommendation'] = getRecommendation(weather['weather_code'])
    return weather

'''Reads the stop info from stops.txt'''
def getStopInfo():
    global stops
    stops = pd.read_csv('stops.txt',index_col='stop_code')

def getRecommendation(code):
    map = {
    200 : 'Get an umbrella!',
    201 : 'Get an umbrella!',
    202 : 'Get an umbrella!',
    210 : 'Get an umbrella!',
    211 : 'Get an umbrella!',
    212 : 'Get an umbrella!',
    221 : 'Get an umbrella!',
    230 : 'Get an umbrella!',
    231 : 'Get an umbrella!',
    232 : 'Get an umbrella!',
    300	: 'Get an umbrella!',
    301	: 'Get an umbrella!',
    302	: 'Get an umbrella!',
    310	: 'Get an umbrella!',
    311	: 'Get an umbrella!',
    312	: 'Get an umbrella!',
    313	: 'Get an umbrella!',
    314	: 'Get an umbrella!',
    321 : 'Get an umbrella!',
    500	: 'Get an umbrella!',
    501	: 'Get an umbrella!',
    502	: 'Get an umbrella!',
    503	: 'Get an umbrella!',
    504	: 'Get an umbrella!',
    511	: 'Get an umbrella!',
    520	: 'Get an umbrella!',
    521	: 'Get an umbrella!',
    522	: 'Get an umbrella!',
    531	: 'Get an umbrella!',
    600	: 'Layer up!',
    601	: 'Layer up!',
    602	: 'Layer up!',
    611	: 'Get an umbrella and layer up!',
    612	: 'Get an umbrella and layer up!',
    615	: 'Get an umbrella and layer up!',
    616	: 'Get an umbrella and layer up!',
    620	: 'Get an umbrella and layer up!',
    621	: 'Get an umbrella and layer up!',
    622	: 'Get an umbrella and layer up!',
    701	: 'Get a light jacket',
    711	: 'Avoid travel',
    721	: 'All set!',
    731	: 'Dusty! Stay safe',
    741	: 'Avoid travel',
    751	: 'All set!',
    761	: 'All set!',
    762	: 'Avoid travel',
    771	: 'Avoid travel',
    781	: 'Avoid travel',
    800	: 'Sunglasses!',
    801	: 'Get an umbrella!',
    802	: 'Get an umbrella!',
    803	: 'Get an umbrella!',
    804	: 'Get an umbrella!',
    }
    return map.get(code)

if __name__ == '__main__':
    app.run()

max_stars_repo_name: vaibhav-walia/roger

max_issues_repo_name: vaibhav-walia/roger

max_forks_repo_name: vaibhav-walia/roger


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

from amical.analysis import pymask

def pymask_grid(input_data, ngrid=40, pa_prior=[0, 360], sep_prior=[0, 100], cr_prior=[1, 150],
                err_scale=1., extra_error_cp=0., ncore=1, verbose=False):
    cpo = pymask.cpo(input_data)
    like_grid = pymask.coarse_grid(cpo, nsep=ngrid, nth=ngrid, ncon=ngrid, thmin=pa_prior[0], thmax=pa_prior[1],
                                   smin=sep_prior[0], smax=sep_prior[1], cmin=cr_prior[0], cmax=cr_prior[1],
                                   threads=ncore, err_scale=err_scale, extra_error=extra_error_cp, verbose=verbose)
    return like_grid

def pymask_mcmc(input_data, initial_guess, niters=1000, pa_prior=[0, 360], sep_prior=None, cr_prior=None,
                err_scale=1, extra_error_cp=0, ncore=1, burn_in=500, walkers=100, display=True,
                verbose=True):
    cpo = pymask.cpo(input_data)
    hammer_data = pymask.hammer(cpo, ivar=initial_guess, niters=niters, model='constant', nwalcps=walkers,
                                sep_prior=sep_prior, pa_prior=pa_prior, crat_prior=cr_prior,
                                err_scale=err_scale, extra_error=extra_error_cp, plot=display,
                                burn_in=burn_in, threads=ncore)

    res_corner = hammer_data[1]
    chain = hammer_data[0]['chain']

    dm = 2.5*np.log10(res_corner['cr'])
    dmm = 2.5*np.log10(res_corner['cr']-res_corner['delcrm'])
    dmp = 2.5*np.log10(res_corner['cr']+res_corner['delcrp'])
    e_dmm = abs(dm - dmm)
    e_dmp = abs(dm - dmp)

    if verbose:
        print('MCMC estimation')
        print('---------------')
        print('Separation = %2.1f +%2.1f/-%2.1f mas' %
              (res_corner['sep'], res_corner['delsepp'], res_corner['delsepm']))
        print('PA = %2.1f +%2.1f/-%2.1f deg' %
              (res_corner['pa'], res_corner['delpap'], res_corner['delpam']))
        print('Contrast Ratio = %2.1f +%2.1f/-%2.1f' %
              (res_corner['cr'], res_corner['delcrp'], res_corner['delcrm']))
        print('dm = %2.2f +%2.2f/-%2.2f mag' % (dm, e_dmp, e_dmm))

    chain_sep = chain[:, :, 0].T
    chain_th = chain[:, :, 1].T
    chain_cr = chain[:, :, 2].T

    if display:
        sep = res_corner['sep']
        pa = res_corner['pa']
        cr = res_corner['cr']
        plt.figure(figsize=(5, 7))
        plt.subplot(3, 1, 1)
        plt.plot(chain_sep, color='grey', alpha=.5)
        plt.plot(len(chain_sep), sep, marker='*', color='#0085ca', zorder=1e3)
        plt.ylabel('Separation [mas]')
        plt.subplot(3, 1, 2)
        plt.plot(chain_th, color='grey', alpha=.5)
        plt.plot(len(chain_sep), pa, marker='*', color='#0085ca', zorder=1e3)
        plt.ylabel('PA [deg]')
        plt.subplot(3, 1, 3)
        plt.plot(chain_cr, color='grey', alpha=.2)
        plt.plot(len(chain_sep), cr, marker='*', color='#0085ca', zorder=1e3)
        plt.xlabel('Step')
        plt.ylabel('CR')
        plt.tight_layout()
        plt.show(block=False)

    res = {'best': {'model': 'binary',
                    'dm': dm,
                    'theta': res_corner['pa'],
                    'sep': res_corner['sep'],
                    'x0': 0,
                    'y0': 0},
           'uncer': {'dm_p': e_dmp,
                     'dm_m': e_dmm,
                     'theta_p': res_corner['delpap'],
                     'theta_m': res_corner['delpam'],
                     'sep_p': res_corner['delsepp'],
                     'sep_m': res_corner['delsepm']},
           }
    return res

def pymask_cr_limit(input_data, nsim=100, err_scale=1, extra_error_cp=0, ncore=1, cmax=500,
                    nsep=60, ncrat=60, nth=30, smin=20, smax=250, cmin=1.0001):
    cpo = pymask.cpo(input_data)
    lims_data = pymask.detec_limits(cpo, threads=ncore, nsim=nsim,
                                    nsep=nsep, ncon=ncrat, nth=nth,
                                    smax=smax, cmax=cmax, cmin=cmin,
                                    smin=smin,
                                    err_scale=err_scale, extra_error=extra_error_cp)

    limits = lims_data['limits']
    seps = lims_data['seps']
    crats = lims_data['cons']
    crat_limits = 0*seps

    # Loop through seps and find the highest contrast ratio that would be detectable
    for sep_ix in range(len(seps)):
        would_detec = limits[:, sep_ix] >= 0.9973
        if np.sum(would_detec) >= 1:
            threesig_lim = np.max(crats[would_detec])
            crat_limits[sep_ix] = threesig_lim  # fivesig_lim
        else:
            crat_limits[sep_ix] = 1.
    con_limits = 2.5*np.log10(crat_limits)

    plt.figure()
    plt.plot(seps, con_limits)
    plt.xlabel('Separation [mas]')
    plt.ylabel('$\Delta \mathrm{Mag}_{3\sigma}$')
    plt.title("PYMASK: flux ratio for 3$\sigma$ detection")
    plt.ylim(plt.ylim()[1], plt.ylim()[0])  # -- rreverse plot
    plt.tight_layout()
    res = {'r': seps,
           'cr_limit': con_limits,
           'lims_data': lims_data
           }
    return res


text:
import numpy as np
import math
from orf import ORF
from utils import dataRange
from firstSliceDataset import dataset_ignit
#import matplotlib.pyplot as plt
from PIL import Image
import cv2
from tqdm import tqdm
from firstSliceSegmentation import firstSliceSegmentation
import datasetManagement
from propagation import propagation

train_data = datasetManagement.getData()

print(train_data)


text:
# cython: language_level=3
# distutils: language = c++
# -*- coding: utf-8 -*-
# *****************************************************************************
# Copyright (c) 2016-2020, Intel Corporation
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# - Redistributions of source code must retain the above copyright notice,
#   this list of conditions and the following disclaimer.
# - Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
# THE POSSIBILITY OF SUCH DAMAGE.
# *****************************************************************************

"""
Interface of the searching function of the dpnp

Notes
-----
This module is a face or public interface file for the library
it contains:
 - Interface functions
 - documentation for the functions
 - The functions parameters check

"""

from dpnp.dpnp_algo import *
from dpnp.dpnp_utils import *

import dpnp
import numpy

__all__ = [
    'argmax',
    'argmin',
    'searchsorted',
    'where'
]

def argmax(x1, axis=None, out=None):
    """
    Returns the indices of the maximum values along an axis.

    For full documentation refer to :obj:`numpy.argmax`.

    Limitations
    -----------
    Input array is supported as :obj:`dpnp.ndarray`.
    Otherwise the function will be executed sequentially on CPU.
    Parameter ``axis`` is supported only with default value ``None``.
    Parameter ``out`` is supported only with default value ``None``.
    Input array data types are limited by supported DPNP :ref:`Data types`.

    See Also
    --------
    :obj:`dpnp.argmin` : Returns the indices of the minimum values along an axis.
    :obj:`dpnp.amax` : The maximum value along a given axis.
    :obj:`dpnp.unravel_index` : Convert a flat index into an index tuple.
    :obj:`dpnp.take_along_axis` : Apply ``np.expand_dims(index_array, axis)``
                                  from argmax to an array as if by calling max.

    Notes
    -----
    In case of multiple occurrences of the maximum values, the indices
    corresponding to the first occurrence are returned.

    Examples
    --------
    >>> import dpnp as np
    >>> a = np.arange(6).reshape((2, 3)) + 10
    >>> a.shape
    (2, 3)
    >>> [i for i in a]
    [10, 11, 12, 13, 14, 15]
    >>> np.argmax(a)
    5

    """

    x1_desc = dpnp.get_dpnp_descriptor(x1)
    if x1_desc:
        if axis is not None:
            pass
        elif out is not None:
            pass
        else:
            result_obj = dpnp_argmax(x1_desc).get_pyobj()
            result = dpnp.convert_single_elem_array_to_scalar(result_obj)

            return result

    return call_origin(numpy.argmax, x1, axis, out)

def argmin(x1, axis=None, out=None):
    """
    Returns the indices of the minimum values along an axis.

    For full documentation refer to :obj:`numpy.argmin`.

    Limitations
    -----------
    Input array is supported as :obj:`dpnp.ndarray`.
    Otherwise the function will be executed sequentially on CPU.
    Parameter ``axis`` is supported only with default value ``None``.
    Parameter ``out`` is supported only with default value ``None``.
    Input array data types are limited by supported DPNP :ref:`Data types`.

    See Also
    --------
    :obj:`dpnp.argmax` : Returns the indices of the maximum values along an axis.
    :obj:`dpnp.amin` : The minimum value along a given axis.
    :obj:`dpnp.unravel_index` : Convert a flat index into an index tuple.
    :obj:`dpnp.take_along_axis` : Apply ``np.expand_dims(index_array, axis)``
                                  from argmin to an array as if by calling min.

    Notes
    -----
    In case of multiple occurrences of the minimum values, the indices
    corresponding to the first occurrence are returned.

    Examples
    --------
    >>> import dpnp as np
    >>> a = np.arange(6).reshape((2, 3)) + 10
    >>> a.shape
    (2, 3)
    >>> [i for i in a]
    [10, 11, 12, 13, 14, 15]
    >>> np.argmin(a)
    0

    """

    x1_desc = dpnp.get_dpnp_descriptor(x1)
    if x1_desc:
        if axis is not None:
            pass
        elif out is not None:
            pass
        else:
            result_obj = dpnp_argmin(x1_desc).get_pyobj()
            result = dpnp.convert_single_elem_array_to_scalar(result_obj)

            return result

    return call_origin(numpy.argmin, x1, axis, out)

def searchsorted(a, v, side='left', sorter=None):
    """
    Find indices where elements should be inserted to maintain order.

    For full documentation refer to :obj:`numpy.searchsorted`.

    """

    return call_origin(numpy.where, a, v, side, sorter)

def where(condition, x=None, y=None):
    """
    Find indices where elements should be inserted to maintain order.

    For full documentation refer to :obj:`numpy.searchsorted`.

    """

    return call_origin(numpy.where, condition, x, y)

max_stars_repo_name: LukichevaPolina/dpnp

max_issues_repo_name: LukichevaPolina/dpnp

max_forks_repo_name: LukichevaPolina/dpnp


text:
#!/usr/bin/env/python3
from scipy.io import wavfile
#import pysptk
import pyreaper
import numpy as np
from pydub import AudioSegment
from pydub.silence import split_on_silence
from pydub.playback import play
import os
import sys

def get_pitch(fs, x):
    pm_times, pm, f0_times, f0, corr = pyreaper.reaper(x, fs)
    values, counts = np.unique(f0, return_counts=True)
    return values[counts.argmax()]

# Set folders and input file here
#file_load = "/home/guy/workspace/wave_detector/dataset4/all_up2.flac"
#output_dir = "/home/guy/workspace/wave_detector/dataset4/all_up"

if len(sys.argv) < 2:
    print("Please provide a file to slice")
    sys.exit(1)
file_load = sys.argv[1]
output_dir = os.path.join(os.path.dirname(file_load), os.path.basename(file_load).split(".")[0])

print("writing file to output:" + output_dir)
os.system("mkdir -p " + output_dir)

#Split audio to bits

def match_target_amplitude(aChunk, target_dBFS):
    ''' Normalize given audio chunk '''
    change_in_dBFS = target_dBFS - aChunk.dBFS
    return aChunk.apply_gain(change_in_dBFS)

song = AudioSegment.from_file(file_load)

#split track where silence is 2 seconds or more and get chunks

chunks = split_on_silence(song,
    # must be silent for at least 0.02 seconds or 20 ms
    min_silence_len=1,

    # consider it silent if quieter than -5 dBFS
    #Adjust this per requirement
    silence_thresh=-30
)

print("chunks: " + str(len(chunks)))

#Process each chunk per requirements
for i, chunk in enumerate(chunks):
    ##Create 0.5 seconds silence chunk
    #silence_chunk = AudioSegment.silent(duration=0)

    #Add  0.5 sec silence to beginning and end of audio chunk
    #audio_chunk = silence_chunk + chunk + silence_chunk
    audio_chunk = chunk

    #Normalize each audio chunk
    normalized_chunk = match_target_amplitude(audio_chunk, -20.0)

    #Export audio chunk with new bitrate
    #print("exporting chunk{0}.mp3".format(i) )
    #normalized_chunk.export(".//chunk{0}.mp3".format(i), bitrate='192k', format="mp3")
    print("exporting chunk{0}.wav".format(i) )

    samples = normalized_chunk.get_array_of_samples()
    #play(normalized_chunk)
    pitch = get_pitch(normalized_chunk.frame_rate, np.array(samples))

    normalized_chunk = normalized_chunk.set_channels(1)
    normalized_chunk.export(os.path.join(output_dir, "chunk{0}-{1}.wav".format(i, str(pitch))), format="wav")

#fs, x = wavfile.read("/home/guy/workspace/wave_detector/audio_analysis/chunk1.wav")

max_stars_repo_name: guysoft/WaveCatchers

max_issues_repo_name: guysoft/WaveCatchers

max_forks_repo_name: guysoft/WaveCatchers

reason: import numpy,from scipy


text:
############################
# Color Range - HSV Unique #
############################
#
#   V1.0
#   29/05/2017
#
#   Implemented by:
#   Yuxi Zhu (matlab) & Thomas Langerak (converted to python)
#   (zhuyuxi1990@gmail.com) & (hello@thomaslangerak.nl)
#
#   Supervisor:
#   Antti Oulasvirta
#
#   This work was funded by Technology Industries of Finland in a three-year
#   project grant on self-optimizing web services. The principal investigator
#   is Antti Oulasvirta of Aalto University (antti.oulasvirta@aalto.fi)
#
###########
# Summary #
###########
#
#   Hasler & Susstrunk validated this metric in their paper. It looks at the number of unique values in the HSV colour
#   space.
#
#############
# Technical #
#############
#
#   Inputs: JPG image (base64)
#   Returns: List of 4 items: Number of Unique HSV (int), Number of Unique Hue (int), Number of Unique Saturation (int), Number of Unique Value (int)
#
##############
# References #
##############
#
#   1.  Hasler, D. and Susstrunk, S. Measuring Colourfuness in Natural Images. (2003).
#
##############
# Change Log #
##############
#
###############
# Bugs/Issues #
###############
#
from skimage import color, util
import collections
import numpy as np
import base64
from PIL import Image
from io import BytesIO

def execute(b64):
    b64 = base64.b64decode(b64)
    b64 = BytesIO(b64)
    img = Image.open(b64)
    img= np.array(img)
    img = util.img_as_ubyte(img)
    img = color.rgb2hsv(img)
    img = img.reshape(-1, 3)
    img = [tuple(l) for l in img]

    hist = collections.Counter(img)
    hist = hist.items()

    hsv_unique = []
    count = []
    h = []
    s = []
    v = []

    for x in range(len(hist)):
        add = [hist[x][0][0], hist[x][0][1], hist[x][0][2]]
        hsv_unique.append(add)
        count.append(hist[x][1])
        h.append(hist[x][0][0])
        s.append(hist[x][0][1])
        v.append(hist[x][0][2])

    # Get all unique values, still has all counts (so no minimal occurence). This probably needs some changing in the future
    h_unique = np.unique(h)
    s_unique = np.unique(s)
    v_unique = np.unique(v)

    new_hsv = []

    # Only often enough occuring values for hsv
    for x in range(len(hsv_unique)):
        if count[x] > 5:
            new_hsv.append(hsv_unique[x])

    result = [len(new_hsv), len(h_unique), len(s_unique), len(v_unique)]

    return result

max_stars_repo_name: kargaranamir/aim

max_issues_repo_name: kargaranamir/aim

max_forks_repo_name: kargaranamir/aim


text:
"""
 Copyright 2018 Johns Hopkins University  (Author: Jesus Villalba)
 Apache 2.0  (http://www.apache.org/licenses/LICENSE-2.0)
"""

import numpy as np

from .roc import compute_rocch, rocch2eer

def compute_dcf(p_miss, p_fa, prior, normalize=True):
    """Computes detection cost function
        DCF = prior*p_miss + (1-prior)*p_fa

    Args:
       p_miss: Vector of miss probabilities.
       p_fa:   Vector of false alarm probabilities.
       prior:  Target prior or vector of target priors.
       normalize: if true, return normalized DCF, else unnormalized.

    Returns:
       Matrix of DCF for each pair of (p_miss, p_fa) and each value of prior.
       [len(prior) x len(p_miss)]
    """

    prior = np.asarray(prior)
    if prior.ndim == 1:
        prior = prior[:, None]

    dcf = prior * p_miss + (1 - prior) * p_fa
    if normalize:
        dcf /= np.minimum(prior, 1 - prior)
    return dcf

def compute_min_dcf(tar, non, prior, normalize=True):
    """Computes minimum DCF
        min_DCF = min_t prior*p_miss(t) + (1-prior)*p_fa(t)
       where t is the decission threshold.

    Args:
      tar: Target scores.
      non: Non-target scores.
      prior: Target prior or vector of target priors.
      normalize: if true, return normalized DCF, else unnormalized.

    Returns:
      Vector Minimum DCF for each prior.
      Vector of P_miss corresponding to each min DCF.
      Vector of P_fa corresponding to each min DCF.
    """

    p_miss, p_fa = compute_rocch(tar, non)
    dcf = compute_dcf(p_miss, p_fa, prior, normalize)
    idx_min_dcf = np.argmin(dcf, axis=-1)
    if dcf.ndim == 1:
        min_dcf = dcf[idx_min_dcf]
        p_miss = p_miss[idx_min_dcf]
        p_fa = p_fa[idx_min_dcf]
    else:
        i1 = np.arange(dcf.shape[0])
        min_dcf = dcf[i1, idx_min_dcf]
        p_miss = p_miss[idx_min_dcf]
        p_fa = p_fa[idx_min_dcf]
    return min_dcf, p_miss, p_fa

def compute_act_dcf(tar, non, prior, normalize=True):
    """Computes actual DCF by making decisions assuming that scores
       are calibrated to act as log-likelihood ratios.

    Args:
      tar: Target scores.
      non: Non-target scores.
      prior: Target prior or vector of target priors.
      normalize: if true, return normalized DCF, else unnormalized.

    Returns:
      Vector actual DCF for each prior.
      Vector of P_miss corresponding to each act DCF.
      Vector of P_fa corresponding to each act DCF.
    """
    prior = np.asarray(prior)

    if prior.ndim == 1:
        assert np.all(
            prior == np.sort(prior, kind="mergesort")
        ), "priors must be in ascending order"
    else:
        prior = prior[None]

    num_priors = len(prior)

    ntar = len(tar)
    nnon = len(non)

    # thresholds
    t = -np.log(prior) + np.log(1 - prior)

    ttar = np.concatenate((t, tar))
    ii = np.argsort(ttar, kind="mergesort")
    r = np.zeros((num_priors + ntar), dtype="int32")
    r[ii] = np.arange(1, num_priors + ntar + 1)
    r = r[:num_priors]
    n_miss = r - np.arange(num_priors, 0, -1)

    tnon = np.concatenate((t, non))
    ii = np.argsort(tnon, kind="mergesort")
    r = np.zeros((num_priors + nnon), dtype="int32")
    r[ii] = np.arange(1, num_priors + nnon + 1)
    r = r[:num_priors]
    n_fa = nnon - r + np.arange(num_priors, 0, -1)

    # n_miss2 = np.zeros((num_priors,), dtype='int32')
    # n_fa2 = np.zeros((num_priors,), dtype='int32')

    # for i in range(len(t)):
    #     n_miss2[i] = np.sum(tar<t[i])
    #     n_fa2[i] = np.sum(non>t[i])

    # assert np.all(n_miss2 == n_miss)
    # assert np.all(n_fa2 == n_fa)
    # print(n_miss)
    # print(n_fa)

    p_miss = n_miss / ntar
    p_fa = n_fa / nnon

    act_dcf = prior * p_miss + (1 - prior) * p_fa
    if normalize:
        act_dcf /= np.minimum(prior, 1 - prior)

    if len(act_dcf) == 1:
        act_dcf = act_dcf[0]

    return act_dcf, p_miss, p_fa

def fast_eval_dcf_eer(tar, non, prior, normalize_dcf=True, return_probs=False):
    """Computes actual DCF, minimum DCF, EER and PRBE all togther

    Args:
      tar: Target scores.
      non: Non-target scores.
      prior: Target prior or vector of target priors.
      normalize_cdf: if true, return normalized DCF, else unnormalized.

    Returns:
      Vector Minimum DCF for each prior.
      Vector Actual DCF for each prior.
      EER value
      PREBP value
    """

    p_miss, p_fa = compute_rocch(tar, non)
    eer = rocch2eer(p_miss, p_fa)

    N_miss = p_miss * len(tar)
    N_fa = p_fa * len(non)
    prbep = rocch2eer(N_miss, N_fa)

    dcf = compute_dcf(p_miss, p_fa, prior, normalize_dcf)
    min_dcf = np.min(dcf, axis=-1)

    act_dcf, act_pmiss, act_pfa = compute_act_dcf(tar, non, prior, normalize_dcf)

    if not return_probs:
        return min_dcf, act_dcf, eer, prbep

    idx = np.argmin(dcf, axis=-1)
    min_pmiss = p_miss[idx]
    min_pfa = p_fa[idx]
    return min_dcf, act_dcf, eer, prbep, min_pmiss, min_pfa, act_pmiss, act_pfa

max_stars_repo_name: hyperion-ml/hyperion

max_issues_repo_name: hyperion-ml/hyperion

max_forks_repo_name: hyperion-ml/hyperion


text:
import numpy as np
import pytest
from .utilities import twenty_one_triangles

# used to parametrize tests for both methods
try:
    import cell_tree2d  # noqa: ignore=F401
    methods = ['simple', 'celltree']
except ImportError:
    # no cell tree -- only test simple
    methods = ['simple']

@pytest.mark.parametrize("method", methods)
def test_single(method, twenty_one_triangles):
    ugrid = twenty_one_triangles
    face = ugrid.locate_faces((4, 6.5), method)
    assert face == 6

@pytest.mark.parametrize("method", methods)
def test_multi(method, twenty_one_triangles):
    ugrid = twenty_one_triangles
    face = ugrid.locate_faces(np.array(((4, 6.5), (7, 2))), method)
    assert (face == np.array((6, 0))).all()

@pytest.mark.parametrize("method", methods)
def test_oob(method, twenty_one_triangles):
    ugrid = twenty_one_triangles
    face = ugrid.locate_faces((0, 0), method)
    assert face == -1
    face = 0
    face = ugrid.locate_faces(np.array(((0, 0),)), method)
    assert np.array_equal(face, np.array((-1, )))

max_stars_repo_name: groutr/gridded

max_issues_repo_name: groutr/gridded

max_forks_repo_name: groutr/gridded