text:
import argparse
import sys
import logging
logger = logging.getLogger(__name__)

from bayes_gain_screens.tomographic_kernel import TomographicKernel
from bayes_gain_screens.utils import make_coord_array
from bayes_gain_screens.plotting import plot_vornoi_map
from bayes_gain_screens.frames import ENU
from h5parm import DataPack
from jaxns.gaussian_process.kernels import RBF
import jax.numpy as jnp
from jax.scipy.linalg import solve_triangular
from jax import jit, random, vmap
from h5parm.utils import create_empty_datapack
import astropy.units as au
import astropy.coordinates as ac
import pylab as plt
import numpy as np

ARRAYS = {'lofar': DataPack.lofar_array_hba}

def get_num_directions(avg_spacing, field_of_view_diameter):
    V = 2.*np.pi*(field_of_view_diameter/2.)**2
    pp = 0.5
    n = -V * np.log(1. - pp) / (avg_spacing/60.)**2 / np.pi / 2.
    n = max(int(n), 50)
    return n

def compute_conditional_moments(kernel:TomographicKernel, X_new, wind_velocity):
    f_K = lambda X1, X2: kernel(X1, X2, bottom=300., width=50., l=4., sigma=1., wind_velocity=wind_velocity)
    K_new_new = f_K(X_new, X_new)
    L_new = jnp.linalg.cholesky(K_new_new + 1e-6*jnp.eye(K_new_new.shape[0]))
    return L_new
    # K_old_old = f_K(X_old, X_old)
    # K_old_new = f_K(X_old, X_new)
    # L = jnp.linalg.cholesky(K_old_old + 1e-6 * jnp.eye(K_old_old.shape[0]))
    # JT = solve_triangular(L, K_old_new, lower=True)
    # C = K_new_new - JT.T @ JT
    # LC = jnp.linalg.cholesky(C + 1e-6*jnp.eye(C.shape[0]))
    # # K_new_old @ (K_old_old)^-1 m(old)
    # # K = L @ L^T
    # # K^-1 = L^-T @ L^-1
    # # (L^-T @ J.T)^T
    # M = solve_triangular(L.T, JT, lower=False)
    # return L_new, LC, M

def main(output_h5parm, ncpu, ra, dec,
         array_name, start_time, time_resolution, duration,
         field_of_view_diameter, avg_direction_spacing, east_wind, north_wind, time_block_size):
    Nd = get_num_directions(avg_direction_spacing, field_of_view_diameter)
    logger.info(f"Number of directions to simulate: {Nd}")
    Nf = 1
    Nt = int(duration / time_resolution) + 1
    time_block_size = min(time_block_size, Nt)
    logger.info(f"Number of times to simulate: {Nt}")
    dp = create_empty_datapack(Nd, Nf, Nt, pols=None,
                          field_of_view_diameter=field_of_view_diameter,
                          start_time=start_time,
                          time_resolution=time_resolution,
                          min_freq=122.,
                          max_freq=166.,
                          array_file=ARRAYS[array_name],
                          phase_tracking=(ra, dec),
                          save_name=output_h5parm,
                          clobber=True)

    with dp:
        dp.current_solset = 'sol000'
        dp.select(pol=slice(0, 1, 1), ant=[0,10,20], time=slice(0,time_block_size))
        axes = dp.axes_tec
        patch_names, directions = dp.get_directions(axes['dir'])
        antenna_labels, antennas = dp.get_antennas(axes['ant'])
        timestamps, times = dp.get_times(axes['time'])
    avg_time = times[len(times) // 2]
    antennas = ac.ITRS(*antennas.cartesian.xyz, obstime=avg_time)
    ref_ant = antennas[0]
    frame = ENU(obstime=avg_time, location=ref_ant.earth_location)
    antennas = antennas.transform_to(frame)
    ref_ant = antennas[0]
    x0 = ref_ant.cartesian.xyz.to(au.km).value
    directions = directions.transform_to(frame)
    t = times.mjd*86400.
    t -= t[0]
    dt = time_resolution
    x = antennas.cartesian.xyz.to(au.km).value.T[1:,:]
    print(x[1]-x[0])
    # x[1,:] = x[0,:]
    # x[1,0] += 0.3
    k = directions.cartesian.xyz.value.T
    logger.info(f"Directions: {directions}")
    logger.info(f"Antennas: {x} {antenna_labels}")
    logger.info(f"Reference Ant: {x0}")
    logger.info(f"Times: {t}")
    Na = x.shape[0]
    logger.info(f"Number of antenna to simulate: {Na}")
    Nd = k.shape[0]
    Nt = t.shape[0]

    #m(X_new) = K(X_new, X_old) @ (K(X_old, X_old))^{-1} m(X_old)
    #K(X_new, X_new) = K(X_new, X_new) - K(X_new, X_old) @ (K(X_old, X_old))^{-1} K(X_old, X_new)

    wind_vector = jnp.asarray([east_wind, north_wind, 0.])/1000.#km/s

    X = make_coord_array(x, k, t[:,None], flat=True)#N,7

    logger.info(f"Sampling {X.shape[0]} new points.")
    kernel = TomographicKernel(x0, x0, RBF(), S_marg=25, compute_tec=True)
    L = jit(compute_conditional_moments, static_argnums=[0])(kernel, X, wind_vector)

    dtec = L @ random.normal(random.PRNGKey(24532), shape=(L.shape[0],1))
    dtec = dtec.reshape((Na, Nd, time_block_size)).transpose((1,0,2))

    with dp:
        dp.select(pol=slice(0, 1, 1), ant=[10, 50], time=slice(0,time_block_size))
        dp.tec = np.asarray(dtec[None, ...])

    fig, axs = plt.subplots(Na,time_block_size, sharex=True, sharey=True, figsize=(2*time_block_size,2*Na))

    for a in range(Na):
        for i in range(time_block_size):
            ax = axs[a][i]
            ax = plot_vornoi_map(k[:, 0:2], dtec[:, a, i], ax=ax, colorbar=False)
            if a == (Na-1):
                ax.set_xlabel(r"$k_{\rm east}$")
            if i == 0:
                ax.set_ylabel(r"$k_{\rm north}$")
            # ax.set_title(f"{} {times[i]}")
            # plt.show()

        # for d in range(Nd):
        #     plt.plot(dtec[d, a,:],alpha=0.3)
        # plt.title(f"{antenna_labels[a]}")
    plt.show()

def debug_main():
    main(output_h5parm='test_datapack.h5',
         ncpu=1,
         ra=120.,
         dec=30.,
         array_name='lofar',
         start_time=None,
         time_resolution=10.,
         duration=600.,
         field_of_view_diameter=4.,
         avg_direction_spacing=8.,
         east_wind=150.,
         north_wind=0.,
         time_block_size=5)

def add_args(parser):
    parser.register("type", "bool", lambda v: v.lower() == "true")
    parser.add_argument('--output_h5parm', help='H5Parm file to file to place the simulated differential TEC',
                        default=None, type=str, required=True)
    parser.add_argument('--ra', help='RA in degrees in ICRS frame.',
                        default=None, type=float, required=True)
    parser.add_argument('--dec', help='DEC in degrees in ICRS frame.',
                        default=None, type=float, required=True)
    parser.add_argument('--array_name', help=f'Name of array, options are {sorted(list(ARRAYS.keys()))}.',
                        default=None, type=float, required=True)
    parser.add_argument('--start_time', help=f'Start time in modified Julian days (mjs/86400).',
                        default=None, type=float, required=True)
    parser.add_argument('--time_resolution', help=f'Temporal resolution in seconds.',
                        default=30., type=float, required=False)
    parser.add_argument('--time_resolution', help=f'Temporal duration in seconds.',
                        default=30., type=float, required=False)
    parser.add_argument('--field_of_view_diameter', help=f'Diameter of field of view in degrees.',
                        default=4., type=float, required=False)
    parser.add_argument('--avg_direction_spacing', help=f'Average spacing between directions in arcmin.',
                        default=6., type=float, required=False)
    parser.add_argument('--east_wind', help=f'Velocity of wind to the east at 100km in m/s.',
                        default=-200., type=float, required=False)
    parser.add_argument('--north_wind', help=f'Velocity of wind to the north at 100km in m/s.',
                        default=0., type=float, required=False)
    parser.add_argument('--ncpu', help='Number of CPUs.',
                        default=1, type=int, required=True)
    parser.add_argument('--time_block_size', help='Number of time steps to simulate at once (must be >= 2).',
                        default=2, type=int, required=True)

if __name__ == '__main__':
    if len(sys.argv) == 1:
        debug_main()
        exit(0)
    parser = argparse.ArgumentParser(
        description='Infers the value of DTEC and a constant over a screen.',
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    add_args(parser)
    flags, unparsed = parser.parse_known_args()
    logger.info("Running with:")
    for option, value in vars(flags).items():
        logger.info("\t{} -> {}".format(option, value))
    main(**vars(flags))

reason: import numpy,import astropy,import jax,from jax


text:
import time

import keras
import matplotlib
import numpy
from matplotlib import pyplot

def plot_learning_curves(loss, val_loss):
    pyplot.plot(numpy.arange(len(loss)) + 0.5, loss, "b.-", label="Training loss")
    pyplot.plot(numpy.arange(len(val_loss)) + 1, val_loss, "r.-", label="Validation loss")
    pyplot.gca().xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
    pyplot.axis([1, 20, 0, 0.05])
    pyplot.legend(fontsize=14)
    pyplot.xlabel("Epochs")
    pyplot.ylabel("Loss")
    pyplot.grid(True)

def plot_series(series, y=None, y_pred=None, x_label="$t$", y_label="$x(t)$", n_steps=50):
    pyplot.plot(series, ".-")
    if y is not None:
        pyplot.plot(n_steps, y, "bx", markersize=10)
    if y_pred is not None:
        pyplot.plot(n_steps, y_pred, "ro")
    pyplot.grid(True)
    if x_label:
        pyplot.xlabel(x_label, fontsize=16)
    if y_label:
        pyplot.ylabel(y_label, fontsize=16, rotation=0)
    pyplot.hlines(0, 0, 100, linewidth=1)
    pyplot.axis([0, n_steps + 1, -1, 1])

def plot_multiple_forecasts(X, Y, Y_pred):
    n_steps = X.shape[1]
    ahead = Y.shape[1]
    plot_series(X[0, :, 0])
    pyplot.plot(numpy.arange(n_steps, n_steps + ahead), Y[0, :, 0], "ro-", label="Actual")
    pyplot.plot(numpy.arange(n_steps, n_steps + ahead), Y_pred[0, :, 0], "bx-", label="Forecast", markersize=10)
    pyplot.axis([0, n_steps + ahead, -1, 1])
    pyplot.legend(fontsize=14)
    pyplot.show()

def generate_time_series(batch_size, n_steps):
    freq1, freq2, offsets1, offsets2 = numpy.random.rand(4, batch_size, 1)
    time_var = numpy.linspace(0, 1, n_steps)
    series = 0.5 * numpy.sin((time_var - offsets1) * (freq1 * 10 + 10))   # wave 1
    series += 0.2 * numpy.sin((time_var - offsets2) * (freq2 * 20 + 20))  # + wave 2
    series += 0.1 * (numpy.random.rand(batch_size, n_steps) - 0.5)    # + noise
    return series[..., numpy.newaxis].astype(numpy.float32)

def create_data_set():
    n_steps = 50
    series = generate_time_series(10000, n_steps + 1)
    X_train, y_train = series[:7000, :n_steps], series[:7000, -1]
    X_valid, y_valid = series[7000:9000, :n_steps], series[7000:9000, -1]
    X_test, y_test = series[9000:, :n_steps], series[9000:, -1]

    return X_train, y_train, X_valid, y_valid, X_test, y_test

def baseline_last_value():
    X_train, y_train, X_valid, y_valid, X_test, y_test = create_data_set()
    y_pred = X_valid[:, -1]
    print(numpy.mean(keras.losses.mean_squared_error(y_valid, y_pred)))

def baseline_linear_regression():
    X_train, y_train, X_valid, y_valid, X_test, y_test = create_data_set()
    model = keras.models.Sequential([
        keras.layers.Flatten(input_shape=[50, 1]),
        keras.layers.Dense(1)
    ])
    model.compile(optimizer='adam', loss='mse')
    model.fit(X_train, y_train, epochs=20)
    y_pred = model.predict(X_valid)
    print(numpy.mean(keras.losses.mean_squared_error(y_valid, y_pred)))

def simple_rnn():
    X_train, y_train, X_valid, y_valid, X_test, y_test = create_data_set()
    model = keras.models.Sequential([
        keras.layers.SimpleRNN(1, input_shape=[None, 1])
    ])
    model.compile(optimizer='adam', loss='mse')
    model.fit(X_train, y_train, epochs=20)
    y_pred = model.predict(X_valid)
    print(model.summary())
    print(numpy.mean(keras.losses.mean_squared_error(y_valid, y_pred)))

def deep_rnn():
    X_train, y_train, X_valid, y_valid, X_test, y_test = create_data_set()
    model = keras.models.Sequential([
        keras.layers.SimpleRNN(20, return_sequences=True, input_shape=[None, 1]),
        keras.layers.SimpleRNN(20, return_sequences=True),
        keras.layers.SimpleRNN(1)
    ])
    model.compile(optimizer='adam', loss='mse')
    print(model.summary())
    # model.fit(X_train, y_train, epochs=20)
    # y_pred = model.predict(X_valid)
    # print(numpy.mean(keras.losses.mean_squared_error(y_valid, y_pred)))

def deep_rnn_dense():
    X_train, y_train, X_valid, y_valid, X_test, y_test = create_data_set()
    model = keras.models.Sequential([
        keras.layers.SimpleRNN(20, return_sequences=True, input_shape=[None, 1]),
        keras.layers.SimpleRNN(20),
        keras.layers.Dense(1)
    ])
    model.compile(optimizer='adam', loss='mse')
    print(model.summary())
    model.fit(X_train, y_train, epochs=20)
    y_pred = model.predict(X_valid)
    print(numpy.mean(keras.losses.mean_squared_error(y_valid, y_pred)))

def deep_rnn_dense_multiple():
    X_train, y_train, X_valid, y_valid, X_test, y_test = create_data_set()
    model = keras.models.Sequential([
        keras.layers.SimpleRNN(20, return_sequences=True, input_shape=[None, 1]),
        keras.layers.SimpleRNN(20),
        keras.layers.Dense(1)
    ])
    model.compile(optimizer='adam', loss='mse')
    model.fit(X_train, y_train, epochs=20)

    n_steps = 50
    series = generate_time_series(1, n_steps + 10)
    X_new, Y_new = series[:, :n_steps], series[:, n_steps:]
    X = X_new
    for step_ahead in range(10):
        y_pred_one = model.predict(X[:, step_ahead:])[:, numpy.newaxis, :]
        X = numpy.concatenate([X, y_pred_one], axis=1)

    Y_pred = X[:, n_steps:]

    plot_multiple_forecasts(X, Y_new, Y_pred)

def main():
    # print(generate_time_series(10, 5))
    # baseline_last_value()
    # baseline_linear_regression()
    # simple_rnn()
    # deep_rnn()
    # deep_rnn_dense()
    # deep_rnn_dense_multiple()
    X_train, y_train, X_valid, y_valid, X_test, y_test = create_data_set()
    print('X_train', X_train.shape)
    print('y_train', y_train.shape)
    print('X_test', X_test.shape)
    print('y_test', y_test.shape)

start = time.time()
main()
end = time.time()
total_time = end - start
print("%s: Total time = %f seconds" % (time.strftime("%Y/%m/%d-%H:%M:%S"), total_time))

ext: py

max_stars_repo_name: franpena-kth/learning-deep-learning

max_issues_repo_name: franpena-kth/learning-deep-learning

max_forks_repo_name: franpena-kth/learning-deep-learning

reason: import numpy


text:
import logging
import scipy as sp

class HammingWindow(object):
    _windows = {}

    @classmethod
    def processBuffer(cls, buf):
        l = len(buf)
        if l not in cls._windows:
            logging.debug("HammingWindow: Generate window for len {}".format(l))
            cls._windows[l] = [  (0.54 - 0.46 * sp.cos(2 * sp.pi * i / (l - 1))) for i in range(l)]

        buf.samples *= cls._windows[l]

reason: import scipy


text:
# Copyright (c) 2011, Roger Lew [see LICENSE.txt]
# This software is funded in part by NIH Grant P20 RR016454.

"""used to estimate RMSD compared to the C-H algorithm for parameter
optimization"""

import math
import time
import numpy as np
from qsturng import qsturng, v_keys
from make_tbls import R

t0 = time.time()
inf = float('inf')

def read_ch(fname):
    with open(fname) as f:
        lines = f.readlines()

    ps,qs = zip(*[L.split(',') for L in lines])

    return map(float, ps), map(float, qs)

ps = np.linspace(.5, .999, 100)
rms = 0.
##for r in [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,
##          22,23,24,25,26,27,28,29,30,35,40,50,60,70,80,90,100,200]:
i=0
nan_count = 0
hmm = 0
error = []
for r in [r for r in [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
          25,30,35,40,50,60,70,80] if r not in R.keys()]:
    for v in [v for v in [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,
              20,22,24,26,30,35,40,50,60,90,120,240,480,inf] if v not in v_keys]:
        ps_ch, qs_ch = read_ch('CH_r=%i,v=%.0f.dat'%(r,v))

        # There are cases where R's qtukey algorithm fails to
        # converge and results in NaN
        for ch, gl, p in zip(qs_ch[:-10],
                             [qsturng(p,r,v) for p in ps[:-10]],
                             ps[:-10]):

##            print ch, gl
            if not math.isnan(ch):

                rms += (ch-gl)**2
                ne = (gl-ch)/ch

                if abs(ne) > .05:
                    print ne,r,v, ch,p

                error.append(ne)

            else:
                nan_count += 1
        i += 1
##        print r,v,rms

rms = math.sqrt(rms / (i*100.))

print rms

print time.time() - t0
import pylab
import numpy as np
from random import shuffle
error = np.array(error)
sq_err = error**2

indices = range(len(error))
sem=[]
for i in range(100):
    shuffle(indices)
    sem.append(np.std(error[indices[:100]])/10.)

pylab.figure()
pylab.title(r'$\mathrm{Bootstrap} \; \mathrm{SEM}$')
pylab.hist(sem)
xticks=pylab.xticks()[0]
pylab.xticks(xticks,[r'$%.1e$'%t for t in xticks])
yticks=pylab.yticks()[0]
pylab.yticks(yticks,[r'$%i$'%t for t in yticks])
pylab.savefig('interp_bootstrap_sd.png')
pylab.close()

mu=[]
for i in range(100):
    shuffle(indices)
    mu.append(math.sqrt(np.sum(sq_err[indices[:100]])/100))

pylab.figure()
pylab.title(r'$\mathrm{Bootstrap} \; \mathrm{RMSD}$')
pylab.hist(sem)
xticks=pylab.xticks()[0]
pylab.xticks(xticks,[r'$%.1e$'%t for t in xticks])
yticks=pylab.yticks()[0]
pylab.yticks(yticks,[r'$%i$'%t for t in yticks])
pylab.savefig('interp_bootstrap_rms.png')
pylab.close()

max_stars_repo_name: rogerlew/qsturng-py

max_issues_repo_name: rogerlew/qsturng-py

max_forks_repo_name: rogerlew/qsturng-py


text:
from pathlib import Path

import cv2
import torch.nn.functional as F
from tqdm import tqdm
import numpy as np

def resize_like(x, target, mode="bilinear"):
    return F.interpolate(x, target.shape[-2:], mode=mode, align_corners=False)

def list2nparray(lst, dtype=None):
    """fast conversion from nested list to ndarray by pre-allocating space"""
    if isinstance(lst, np.ndarray):
        return lst
    assert isinstance(lst, (list, tuple)), "bad type: {}".format(type(lst))
    assert lst, "attempt to convert empty list to np array"
    if isinstance(lst[0], np.ndarray):
        dim1 = lst[0].shape
        assert all(i.shape == dim1 for i in lst)
        if dtype is None:
            dtype = lst[0].dtype
            assert all(i.dtype == dtype for i in lst), "bad dtype: {} {}".format(
                dtype, set(i.dtype for i in lst)
            )
    elif isinstance(lst[0], (int, float, complex, np.number)):
        return np.array(lst, dtype=dtype)
    else:
        dim1 = list2nparray(lst[0])
        if dtype is None:
            dtype = dim1.dtype
        dim1 = dim1.shape
    shape = [len(lst)] + list(dim1)
    rst = np.empty(shape, dtype=dtype)
    for idx, i in enumerate(lst):
        rst[idx] = i
    return rst

def get_img_list(path):
    if Path(path).is_file():
        return [Path(path)]
    else:
        return (
            sorted(list(Path(path).glob("*.png")))
            + sorted(list(Path(path).glob("*.jpg")))
            + sorted(list(Path(path).glob("*.jpeg")))
        )

def gen_miss(img, mask, output):

    imgs = get_img_list(img)
    masks = get_img_list(mask)
    print("Total images:", len(imgs), len(masks))

    out = Path(output)
    out.mkdir(parents=True, exist_ok=True)

    for i, (img, mask) in tqdm(enumerate(zip(imgs, masks))):
        path = out.joinpath("miss_%04d.png" % (i + 1))
        img = cv2.imread(str(img), cv2.IMREAD_COLOR)
        mask = cv2.imread(str(mask), cv2.IMREAD_GRAYSCALE)
        mask = cv2.resize(mask, img.shape[:2][::-1])
        mask = mask[..., np.newaxis]
        miss = img * (mask > 127) + 255 * (mask <= 127)
        cv2.imwrite(str(path), miss)

def merge_imgs(dirs, output, row=1, gap=2, res=512):

    image_list = [get_img_list(path) for path in dirs]
    img_count = [len(image) for image in image_list]
    print("Total images:", img_count)
    assert min(img_count) > 0, "Please check the path of empty folder."

    output_dir = Path(output)
    output_dir.mkdir(parents=True, exist_ok=True)

    n_img = len(dirs)
    row = row
    column = (n_img - 1) // row + 1
    print("Row:", row)
    print("Column:", column)

    for i, unit in tqdm(enumerate(zip(*image_list))):
        name = output_dir.joinpath("merge_%04d.png" % i)
        merge = (
            np.ones(
                [res * row + (row + 1) * gap, res * column + (column + 1) * gap, 3],
                np.uint8,
            )
            * 255
        )
        for j, img in enumerate(unit):
            r = j // column
            c = j - r * column
            img = cv2.imread(str(img), cv2.IMREAD_COLOR)
            if img.shape[:2] != (res, res):
                img = cv2.resize(img, (res, res))
            start_h, start_w = (r + 1) * gap + r * res, (c + 1) * gap + c * res
            merge[start_h : start_h + res, start_w : start_w + res] = img
        cv2.imwrite(str(name), merge)


text:
'''
Miscellaneous data generator utilities.

Copyright (C) 2018 Pierluigi Ferrari

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

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

from __future__ import division
import numpy as np

def apply_inverse_transforms(y_pred_decoded, inverse_transforms):
    '''
    将变换后的图像的预测值，转换回原始图像对应的预测值。

    Arguments:
        y_pred_decoded (list or array):(num_predictions, 6)/(batch_size, num_predictions, 6)`.
            最后一维 类别ID, 类别分数, 四个坐标
        inverse_predictions (list): A nested list of length `batch_size` that contains
            for each batch item a list of functions that take one argument (one element
            of `y_pred_decoded` if it is a list or one slice along the first axis of
            `y_pred_decoded` if it is an array) and return an output of the same shape
            and data type.

    Returns:
        The transformed predictions, which have the same structure as `y_pred_decoded`.
    '''

    if isinstance(y_pred_decoded, list):

        y_pred_decoded_inv = []

        for i in range(len(y_pred_decoded)):
            y_pred_decoded_inv.append(np.copy(y_pred_decoded[i]))
            if y_pred_decoded_inv[i].size > 0: # If there are any predictions for this batch item.
                for inverter in inverse_transforms[i]:
                    if not (inverter is None):
                        y_pred_decoded_inv[i] = inverter(y_pred_decoded_inv[i])

    elif isinstance(y_pred_decoded, np.ndarray):

        y_pred_decoded_inv = np.copy(y_pred_decoded)

        for i in range(len(y_pred_decoded)):
            if y_pred_decoded_inv[i].size > 0: # If there are any predictions for this batch item.
                for inverter in inverse_transforms[i]:
                    if not (inverter is None):
                        y_pred_decoded_inv[i] = inverter(y_pred_decoded_inv[i])

    else:
        raise ValueError("`y_pred_decoded` must be either a list or a Numpy array.")

    return y_pred_decoded_inv

max_stars_repo_name: CVBase-Bupt/SSD-keras

max_issues_repo_name: CVBase-Bupt/SSD-keras

max_forks_repo_name: CVBase-Bupt/SSD-keras


text:
from __future__ import print_function, absolute_import
from six import iteritems, itervalues, iterkeys
from six.moves import range

import tables
import numpy as np

from card_table import CardTable, TableDef
from data_helper import DataHelper

class Parameter(object):
    def __init__(self, h5n, input):
        self._h5n = h5n
        self._input = input

        self.nlparm = NLPARM(self._h5n, self)

    def path(self):
        return self._input.path() + ['PARAMETER']

    def read(self):
        for key, item in iteritems(self.__dict__):
            if key.startswith('_'):
                continue
            try:
                item.read()
            except AttributeError:
                pass

########################################################################################################################

class NLPARM(CardTable):
    table_def = TableDef.create('/NASTRAN/INPUT/PARAMETER/NLPARM')

    @classmethod
    def from_bdf(cls, cards):
        card_ids = sorted(cards.keys())

        data = np.empty(len(card_ids), dtype=cls.table_def.dtype)

        sid = data['SID']
        ninc = data['NINC']
        dt = data['DT']
        kmethod = data['KMETHOD']
        kstep = data['KSTEP']
        maxiter = data['MAXITER']
        conv = data['CONV']
        intout = data['INTOUT']
        epsu = data['EPSU']
        epsp = data['EPSP']
        epsw = data['EPSW']
        maxdiv = data['MAXDIV']
        maxqn = data['MAXQN']
        maxls = data['MAXLS']
        fstress = data['FSTRESS']
        lstol = data['LSTOL']
        maxbis = data['MAXBIS']
        maxr = data['MAXR']
        rtolb = data['RTOLB']
        miniter = data['MINITER']

        # TODO: check that kmethod is correct
        _kmethod = {'AUTO': 0, 'ITER': 1, 'SEMI': 2, 'FNT': 3, 'PFNT': 4, None: 0, '': 0}

        # TODO: check that conv is correct; what are the allowed combinations?
        _conv = {}

        # TODO: check that intout is correct
        _intout = {'YES': 1, 'NO': 0, 'ALL': 3, None: DataHelper.default_int, '': DataHelper.default_int}

        i = -1
        for card_id in card_ids:
            i += 1
            card = cards[card_id]

            sid[i] = card.nlparm_id
            ninc[i] = card.ninc
            dt[i] = card.dt
            kmethod[i] = _kmethod[card.kmethod]
            kstep[i] = card.kstep
            maxiter[i] = card.max_iter
            conv[i] = DataHelper.unknown_int
            if isinstance(card.int_out, int):
                intout[i] = card.int_out
            else:
                intout[i] = _intout[card.int_out]
            epsu[i] = card.eps_u
            epsp[i] = card.eps_p
            epsw[i] = card.eps_w
            maxdiv[i] = card.max_div
            maxqn[i] = card.max_qn
            maxls[i] = card.max_ls
            fstress[i] = card.fstress
            lstol[i] = card.ls_tol
            maxbis[i] = card.max_bisect
            maxr[i] = card.max_r
            rtolb[i] = card.rtol_b
            try:
                miniter[i] = card.min_iter
            except AttributeError:
                miniter[i] = DataHelper.default_int

        result = {'IDENTITY': data}

        return result


text:
from typing import Union
import numpy as np
import pandas as pd

from ..utils import np_sample

def sample_data(X: Union[pd.DataFrame, np.ndarray],
                n=None,
                frac=None,
                random_state=None,
                replace=False,
                **kwargs) -> Union[pd.DataFrame, np.ndarray]:
    assert frac is None or n is None, "Cannot specify both `n` and `frac`"
    assert not (frac is None and n is None), "One of `n` or `frac` must be passed."

    input_is_pd = isinstance(X, pd.DataFrame)

    if not input_is_pd:
        X = pd.DataFrame(X)

    sample = X.sample(
        n=n,
        frac=frac,
        random_state=random_state,
        replace=replace,
        **kwargs,
    )

    if not input_is_pd:
        return sample.values
    else:
        return sample

max_stars_repo_name: emanuele-albini/emutils

max_issues_repo_name: emanuele-albini/emutils

max_forks_repo_name: emanuele-albini/emutils


text:
"""
test_classification.py tests the success of my classification
method Bubble.classify() by using a modified version added to
the test version of the object, BubbleTest.classify_test(), which
indicates which level of classification the object passed (looks like
a bubble) and failed (does not look like a bubble).

***Run this on many datasets until you hone in the classification parameters***
OR, write an ML algo to perform this (requires efficient labeling pipeline)

Author: Andy Ylitalo
Date: July 5, 2021
"""

import os
import pickle as pkl
import numpy as np
import glob

import sys
sys.path.append('../src/')
import config as cfg
import cvimproc.basic as basic
import genl.conversions as conv
import genl.readin as readin
import genl.main_helper as mh

def classify_test(obj, max_aspect_ratio=10, min_solidity=0.85,
            min_size=50, max_orientation=np.pi/8, circle_aspect_ratio=1.1,
            n_consec=3):
    """
    Classifies Bubble objects into different categories. Still testing to see
    what classifications are most important.

    Parameters are currently guesses. ML would be a better method for selecting them.

    Turn this into ML with the following labels (sequences of):
    - Aspect Ratio
    - Orientation
    - Area
    - Rate of area growth
    - Speed along flow direction
    - Solidity
    - *averages thereof
    - predicted max speed
    - width of inner stream
    - height

    ***Ask Chris about this! TODO
    ***CHALLENGE: how do I classify large, elongated bubbles separately? TODO
    """
    # creates dictionary of classifications
    classifications = {'solid' : None, 'circular' : None, 'oriented' : None,
                        'consecutive' : None, 'inner stream' : None,
                        'exited' : None, 'growing' : None}

    # computes average speed [m/s] if not done so already
    if obj.get_props('average flow speed [m/s]') == None:
        obj.process_props()

    # determines which frames are worth looking at (large enough object, not on border)
    worth_checking = np.logical_and(
                            np.asarray(obj.get_props('area')) > min_size,
                            np.logical_not(np.asarray(obj.get_props('on border')))
                            )

    # first checks if object is an artifact (marked by inner_stream = -1)
    # not convex enough?
    # TODO -- scale solidity requirement by size (currently lets through many non-solid large bubbles)
    classifications['solid'] = np.logical_not(np.logical_and(
                np.asarray(obj.get_props('solidity')) < min_solidity,
                worth_checking ))

    # too oblong?
    classifications['circular'] = np.logical_not(np.logical_and(
                np.asarray(obj.get_props('aspect ratio')) > max_aspect_ratio,
                np.asarray([row_max - row_min for row_min, _, row_max, _ in obj.get_props('bbox')]) \
                < obj.get_metadata('R_i')*conv.m_2_um*obj.get_metadata('pix_per_um') ))

    # oriented off axis (orientation of skimage.regionprops is CCW from up direction)
    classifications['oriented'] = np.logical_not(np.logical_and(
                np.logical_and(
                    np.abs(np.asarray(obj.get_props('orientation')) + np.pi/2) > max_orientation,
                    np.abs(np.asarray(obj.get_props('orientation')) - np.pi/2) > max_orientation),
                worth_checking ))

    # checks that the bubble appears in a minimum number of consecutive frames
    if obj.appears_in_consec_frames(n_consec=n_consec):
        classifications['consecutive'] = np.ones([len(obj.get_props('frame'))])
    else:
        classifications['consecutive'] = np.zeros([len(obj.get_props('frame'))])

    # flowing backwards
    classifications['inner stream'] = np.asarray(obj.get_props('flow speed [m/s]')) > obj.get_metadata('v_interf')

    # exits properly from downstream side of frame
    # farthest right column TODO -- determine coordinate of bbox from flow dir
    col_max = obj.get_props('bbox')[-1][-1]
    # average flow speed [m/s]
    average_flow_speed_m_s = obj.get_props('average flow speed [m/s]')
    # number of pixels traveled per frame along flow direction
    if average_flow_speed_m_s is not None:
        pix_per_frame_along_flow = average_flow_speed_m_s * conv.m_2_um * \
                        obj.get_metadata('pix_per_um') / obj.get_metadata('fps')
        # checks if the next frame of the object would reach the border or if the
        # last view of object is already on the border
        if (col_max + pix_per_frame_along_flow < obj.get_metadata('frame_dim')[1]) and \
                not obj.get_props('on border')[-1]:
            classifications['exited'] = np.zeros([len(obj.get_props('frame'))])
        else:
            classifications['exited'] = np.ones([len(obj.get_props('frame'))])
    else:
        classifications['exited'] = np.zeros([len(obj.get_props('frame'))])

    # growing over time
    frames = np.asarray(obj.get_props('frame'))
    areas = np.asarray(obj.get_props('area'))
    if len(frames[worth_checking]) > 1:
        growing = np.polyfit(frames[worth_checking], areas[worth_checking], 1)[0] >= 0
    else:
        growing = True
    if growing:
        classifications['growing'] = np.ones([len(obj.get_props('frame'))])
    else:
        classifications['growing'] = np.zeros([len(obj.get_props('frame'))])

    return classifications

def load_objs(input_filepath, ext, save_dir):
    """
    Loads objects pointed to by input filepath.

    Parameters
    ----------
    input_filepath : string
        Path to input file (stored in `data` subfolder)
    ext : string
        extension of files saved before, without '.' (e.g., 'jpg')

    Returns
    -------
    objs : dictionary
        Dictionary of TrackedObj objects
    """
    # loads parameters
    p = readin.load_params(input_filepath)
    # determines appropriate directories based on parameters
    vid_path, data_path, vid_dir, \
    data_dir, figs_dir, _ = mh.get_paths(p)

    # determines save paths
    p['vid_subdir'] = save_dir
    _, _, _, _, figs_dir, _ = mh.get_paths(p)

    # if saving data, removes existing data
    files_to_remove = glob.glob(os.path.join(figs_dir, '*.{0:s}'.format(ext)))
    for filepath in files_to_remove:
        os.remove(filepath)

    # loads data
    with open(data_path, 'rb') as f:
        data = pkl.load(f)
        objs = data['objects']

    return objs

def superpose_images(obj, classifications, skip_overlaps=False, prop='', color=(255,0,0)):
    """
    Superposes images of an object onto one frame.

    Parameters
    ----------
    obj : TrackedObject
        Object that has been tracked. Must have 'image', 'local centroid',
        'frame_dim', 'bbox', and 'centroid' parameters.

    Returns
    -------
    im : (M x N) numpy array of uint8
        Image of object over time; each snapshot is superposed
        (likely black-and-white)
    """
    # gets dimensions of frame
    M, N = obj.get_metadata('frame_dim')
    # initializes image black
    im = np.zeros([M, N, 3], dtype='uint8')

    # initializes previous bounding box
    bbox_prev = (0,0,0,0)

    # superposes image from each frame
    frame_list = obj.get_props('frame')
    for i, f in enumerate(frame_list):
        # loads bounding box and image within it
        bbox = obj.get_prop('bbox', f)

        # skips images that overlap if requested
        if skip_overlaps:
            if basic.is_overlapping(bbox_prev, bbox):
                continue
            else:
                bbox_prev = bbox

        # loads image and colors it
        im_obj_bw = obj.get_prop('image', f)
        m, n = im_obj_bw.shape
        im_obj = np.zeros([m, n, 3], dtype='uint8')

        # changes color to white if property required
        color_curr = color
        if len(prop) > 0:
            if classifications[prop][i]:
                color_curr = cfg.white
        else:
            color_curr = cfg.white
        # colors image
        im_obj[im_obj_bw, :] = color_curr

        # superposes object image on overall image
        row_min, col_min, row_max, col_max = bbox
        im[row_min:row_max, col_min:col_max, :] = im_obj

    return im

###############################################################################

# OR, DO I WANT TO PROCESS THESE SPECIALLY?

def main():

    # input_subpaths = ['sd301_co2/20210331_45bar/sd301_co2_40000_001_050_0150_95_04_9/testclassify/data/input.txt',
    #                 'sd301_co2/20210207_88bar/sd301_co2_15000_001_100_0335_79_04_10/testclassify/data/input.txt']

    input_subpaths = ['ppg_co2/20210720_70bar/ppg_co2_40000_001-1_050_0228_74_04_11/nofill/data/input.txt']
    ext = 'jpg'
    skip_overlaps = True
    save = True
    prop_list = ['solid', 'circular', 'oriented', ''] # default, no properties
    cmap = {'solid' : (255,0,0), 'circular' : (0,255,0), 'oriented' : (0,0,255),
                        'consecutive' : (255,255,0), 'inner stream' : (255,0,255),
                        'exited' : (0,255,255), 'growing' : (100, 250, 100), '' : (255,255,255)}
    save_dir = 'test/classification'

    input_filepaths = [os.path.join(cfg.output_dir, input_subpath) for input_subpath in input_subpaths]

    # loops through datasets
    for input_filepath in input_filepaths:
        print(input_filepath)

        objs = load_objs(input_filepath, ext, save_dir)

        # loops through objects
        for ID, obj in objs.items():
            # perform test classification
            classifications = classify_test(obj)
            # prints out classifications
            for prop, val in classifications.items():
                if any(np.logical_not(val)):
                    print('Object {0:d} is not {1:s}'.format(ID, prop))

            if save:
                for prop in prop_list:
                    # creates image of bubbles superposed on frame
                    im = superpose_images(obj, classifications, skip_overlaps=skip_overlaps, prop=prop, color=cmap[prop])
                    # saves image
                    basic.save_image(im, os.path.join(cfg.output_dir, save_dir, '{0:s}{1:d}.{2:s}'.format(prop, ID, ext)))

    return

if __name__ == '__main__':
    main()


text:
""" Input/Output helpers

Classes:

* :py:class:`DataReader`: Read a number of imaging and ephys file formats

Functions:

* :py:func:`find_ca_data`: Find trace data under a directory
* :py:func:`save_outfile`: Write out the individual processed traces
* :py:func:`save_final_stats`: Write out the summary stats for all the traces

"""

# Imports
import re
import inspect
import pathlib
from typing import Dict, Callable, Generator, Tuple

# 3rd party
import numpy as np

import pandas as pd

# Our own imports
from .consts import DATA_TYPE

# Column Name Parsers

TIME_KEY1 = re.compile(r'''^
    time\:\:.*time\!\!r
$''', re.IGNORECASE | re.VERBOSE)
AREA_KEY1 = re.compile(r'''^
    channel[0-9]+_R[0-9]+_area\:\:.*area\!\!r
$''', re.IGNORECASE | re.VERBOSE)
MEAN_KEY1 = re.compile(r'''^
    channel[0-9]+_R[0-9]+_intensitymean\:\:.*intensitymean\!\!r
$''', re.IGNORECASE | re.VERBOSE)

TIME_KEYS = [
    TIME_KEY1,
    re.compile(r'^time$', re.IGNORECASE),
    re.compile(r'^.*t\:\:time\!\!.*$', re.IGNORECASE),
]
AREA_KEYS = [
    AREA_KEY1,
    re.compile(r'^area\s*[0-9]+$', re.IGNORECASE),
    re.compile(r'^.*area\:\:area\!\!.*$', re.IGNORECASE),
]
MEAN_KEYS = [
    MEAN_KEY1,
    re.compile(r'^mean\s*[0-9]+$', re.IGNORECASE),
    re.compile(r'^.*value\:\:intensity.*$', re.IGNORECASE),
]
NAME_KEYS = [
    re.compile(r'^.*name\:\:name.*$', re.IGNORECASE),
]

# File Name Patterns

READER_NAME = re.compile(r'^read_infile_(?P<name>[a-z0-9]+)_data$', re.IGNORECASE)

SUFFIXES = ('.csv', '.tsv', '.xls', '.xlsx')

# Classes

class DataReader(object):
    """ Read data from different data types

    Readers for different file types:

    * :py:meth:`read_infile_ca_data`: Data from Zen GCAMP imaging
    * :py:meth:`read_infile_ephys_data`: Data from electrophysiology
    * :py:meth:`read_infile_measure_data`: Data from the Zen measure tool

    Define a new reader:

    .. code-block:: python

        def read_infile_mytype_data(self, infile, **kwargs):
            # read the file here

    To use call:

    .. code-block:: bash

        $ analyze_ca_data.py --data-type mytype /path/to/data

    """

    def __init__(self, data_type: str = DATA_TYPE):
        self.data_type = data_type

    @classmethod
    def get_readers(cls) -> Dict[str, Callable]:
        """ Get the supported readers for this class

        :returns:
            A dictionary of reader_name: reader_method pairs
        """
        readers = {}
        for meth_name, meth in inspect.getmembers(cls(), inspect.ismethod):
            match = READER_NAME.match(meth_name)
            if not match:
                continue
            readers[match.group('name')] = meth
        return readers

    def read_data_frame(self, infile: pathlib.Path, **kwargs) -> pd.DataFrame:
        """ Read in a data frame from a spreadsheet

        :param Path infile:
            The input file to read
        :param \\*\\*kwargs:
            Arguments to pass to the various ``pandas.read_*`` methods
        :returns:
            A DataFrame containing the contents of the file
        """
        infile = pathlib.Path(infile)
        if not infile.is_file():
            raise OSError(f'Cannot find file "{infile}"')

        # Read in the file based on file suffix
        if infile.name.endswith('.csv'):
            return pd.read_csv(str(infile), sep=',', **kwargs)
        elif infile.name.endswith('.tsv'):
            return pd.read_csv(str(infile), sep='\t', **kwargs)
        elif infile.name.endswith(('.xls', '.xlsx')):
            return pd.read_excel(str(infile), **kwargs)
        else:
            raise KeyError(f'No DataFrame reader for file type "{infile}"')

    def read_infile(self, infile: pathlib.Path, **kwargs) -> Tuple[np.ndarray]:
        """ Dispatch based on type

        :param Path infile:
            The input file to read
        :param \\*\\* kwargs:
            The arguments to the individual readers
        """
        flip_signal = kwargs.pop('flip_signal', False)
        level_shift = kwargs.pop('level_shift', None)

        # Tricky method lookup dispatch
        # Call a method named ``read_infile_{data_type}`` and return the results
        time, means = self.get_readers()[self.data_type](infile, **kwargs)

        # Calculate the min/max of each signal
        min_means = np.nanmin(means, axis=0)
        max_means = np.nanmax(means, axis=0)
        rng_means = (max_means - min_means)

        if flip_signal:
            # Flip the signal without level shifting
            means = (means - min_means) / rng_means
            means = (1.0 - means) * rng_means + min_means

        # Level shift all the means to prevent stupid things happening around 0
        level_min = 1.0
        if level_shift in (None, 'min'):
            shift = np.abs(np.nanmin(means, axis=0))
        elif level_shift == 'range':
            shift = np.abs(rng_means.copy())
        elif level_shift == 'one':
            shift = np.ones_like(rng_means)
        elif level_shift == 'zero':
            shift = np.zeros_like(rng_means)
            level_min = 0.0  # hope you know what you're doing...
        else:
            raise ValueError(f'Unknown level shift setting: {level_shift}')

        # After this point, the signal should be everywhere positive and >= level_min
        shift[shift < level_min] = level_min
        means = means - np.nanmin(means, axis=0) + shift

        print(f'Got {means.shape[1]} traces with {means.shape[0]} samples each')
        assert means.shape[0] == time.shape[0]
        return time, means

    # File type readers

    def read_infile_ephys_data(self, infile: pathlib.Path) -> Tuple[np.ndarray]:
        """ Read Mike's data

        Ephys data is a CSV file with or without a header with columns:

        .. code-block:: text

            Time,Mean1,Mean2,Mean3,...
            1,0.1,0.2,0.3,...
            2,0.2,0.3,0.4,...

        All means are assumed aligned in time with no filtering for invalid values.

        Most of these types of files only have one mean/file, which makes it easier.

        :param Path infile:
            The input '.csv' files to read
        :returns:
            (time, means) for each roi in the file
        """
        # Try the first time without a header
        df = self.read_data_frame(infile, header=None)
        time = df.iloc[:, 0]
        if time.dtype not in (np.dtype('float64'), np.dtype('int64')):
            df = self.read_data_frame(infile, header=None, skiprows=1)
        time = df.iloc[:, 0]
        if time.dtype not in (np.dtype('float64'), np.dtype('int64')):
            raise ValueError(f'Cannot understand file header: {infile}')

        # The first vector is time, all subsequent vectors are traces
        time = df.iloc[:, 0].values
        means = df.iloc[:, 1:].values

        assert time.shape[0] == means.shape[0]
        return time, means

    def read_infile_ca_data(self,
                            infile: pathlib.Path,
                            min_samples: int = 10) -> Tuple[np.ndarray]:
        """ Read the Ca data input file

        CA GCAMP imaging has a crazy file structure:

        .. code-block:: text

            WeirdTimeColumn,Area1,Mean1,Total1,Max1,Area2,Mean2,...
            "MS",stupid,number,of,commas,...
            1,57,0.1,inf,inf,58,0.2,...

        If the means are not aligned in time, the time or area values are negative,
        this reader masks those values out as NaNs.

        This loader also tries to filter out super short reads, because they seem to
        be errors too...

        :param Path infile:
            The input '.csv' file to read
        :param int min_samples:
            Minimum number of samples to allow in a time series
        :returns:
            (time, means) for each roi in the file
        """

        # Use a regex on the headers to match time/signal/area
        time_indices = []
        area_indices = []
        mean_indices = []

        df = self.read_data_frame(infile)
        df = df.iloc[1:, :]

        column_matchers = [
            (TIME_KEYS, time_indices),
            (AREA_KEYS, area_indices),
            (MEAN_KEYS, mean_indices),
        ]

        for colidx, column in enumerate(df.columns):
            # Work out column indices by matching on the name
            for colkeys, colindices in column_matchers:
                had_match = False
                for colkey in colkeys:
                    if colkey.match(str(column).strip()):
                        colindices.append(colidx)
                        had_match = True
                        break
                if had_match:
                    break

        # Sanity checks for badly formatted files

        # Time index needs to be unique
        if len(time_indices) == 0:
            err = f'Invalid stat file. Cannot find the time index in: {infile}'
            raise OSError(err)
        elif len(time_indices) > 1:
            err = f'Invalid stat file. Columns {time_indices} all appear to be time records in: {infile}'
            raise OSError(err)
        time_idx = time_indices[0]

        # Area index and mean need to be non-zero length
        if len(area_indices) == 0:
            err = f'Invalid stat file. Cannot find any area columns in: {infile}'
            raise OSError(err)
        if len(mean_indices) == 0:
            err = f'Invalid stat file. Cannot find any signal value columns in: {infile}'
            raise OSError(err)

        # Need the same number of areas as means
        if len(area_indices) != len(mean_indices):
            err = 'Invalid stat file. Got {} unique ROIs but {} signal values in: {}'
            err = err.format(len(area_indices), len(mean_indices), infile)
            raise OSError(err)

        # Make sure we're extracting the columns properly
        # print('Time column:  {}'.format(time_idx))
        # print('Area columns: {}'.format(area_indices))
        # print('Mean columns: {}'.format(mean_indices))

        time = df.iloc[:, time_idx].values.astype(np.float)
        areas = df.iloc[:, area_indices].values
        means = df.iloc[:, mean_indices].values

        # Area == 0 means that the ROI isn't valid for this timepoint
        means[areas <= 0] = np.nan

        # Time < 0 means that this sample isn't valid for this timepoint
        means[time < 0] = np.nan

        # Throw out any columns with too many empty samples
        keep_cols = np.sum(~np.isnan(means), axis=0) >= min_samples
        means = means[:, keep_cols]
        if means.shape[1] < 1:
            err = f'Empty stat file. No ROIs with at least {min_samples} samples in {infile}'
            raise OSError(err)

        return time, means

    def read_infile_measure_data(self,
                                 infile: pathlib.Path,
                                 min_samples: int = 10) -> Tuple[np.ndarray]:
        """ Read the Zen Measure Tool Output file

        CA GCAMP imaging has a crazy file structure:

        .. code-block:: text

            WeirdNameColumn,WeirdTime,WeirdArea,WeirdDiameter
            "MS",stupid,number,of,commas,...
            name1,time1,area1,diameter1,
            name2,time1,area1,diameter1,
            name1,time2,area2,diameter2,
            name2,time2,area2,diameter2,

        If the means are not aligned in time, the time or area values are negative,
        this reader masks those values out as NaNs.

        This loader also tries to filter out super short reads, because they seem to
        be errors too...

        :param Path infile:
            The input '.csv' file to read
        :param int min_samples:
            Minimum number of samples to allow in a time series
        :returns:
            (time, means) for each roi in the file
        """

        # Use a regex on the headers to match time/signal/area
        name_indices = []
        time_indices = []
        area_indices = []
        mean_indices = []

        df = self.read_data_frame(infile)
        df = df.iloc[1:, :]

        column_matchers = [
            (NAME_KEYS, name_indices),
            (TIME_KEYS, time_indices),
            (AREA_KEYS, area_indices),
            (MEAN_KEYS, mean_indices),
        ]

        for colidx, column in enumerate(df.columns):
            # Work out column indices by matching on the name
            for colkeys, colindices in column_matchers:
                had_match = False
                for colkey in colkeys:
                    if colkey.match(str(column).strip()):
                        colindices.append(colidx)
                        had_match = True
                        break
                if had_match:
                    break

        # Sanity checks for badly formatted files

        # For this format, all columns need to be unique
        if len(time_indices) == 0:
            err = f'Invalid stat file. Cannot find the time index in: {infile}'
            raise OSError(err)
        elif len(time_indices) > 1:
            err = f'Invalid stat file. Columns {time_indices} all appear to be time records in: {infile}'
            raise OSError(err)
        time_idx = time_indices[0]

        if len(area_indices) == 0:
            err = f'Invalid stat file. Cannot find any area columns in: {infile}'
            raise OSError(err)
        elif len(area_indices) > 1:
            err = f'Invalid stat file. Columns {area_indices} all appear to be area records in: {infile}'
            raise OSError(err)
        area_idx = area_indices[0]

        if len(mean_indices) == 0:
            err = f'Invalid stat file. Cannot find any signal value columns in: {infile}'
            raise OSError(err)
        elif len(mean_indices) > 1:
            err = f'Invalid stat file. Columns {mean_indices} all appear to be signal value records in: {infile}'
            raise OSError(err)
        mean_idx = mean_indices[0]

        if len(name_indices) == 0:
            err = f'Invalid stat file. Cannot find any name columns in: {infile}'
            raise OSError(err)
        elif len(name_indices) > 1:
            err = f'Invalid stat file. Columns {name_indices} all appear to be name records in: {infile}'
            raise OSError(err)
        name_idx = name_indices[0]

        # Make sure we're extracting the columns properly
        # print('Time column:  {}'.format(time_idx))
        # print('Name column: {}'.format(name_idx))
        # print('Area column: {}'.format(area_idx))
        # print('Mean column: {}'.format(mean_idx))

        # Now restack everything so there is one time column, many signal columns
        names = df.iloc[:, name_idx].values
        times = df.iloc[:, time_idx].values.astype(np.float)
        areas = df.iloc[:, area_idx].values.astype(np.float)
        means = df.iloc[:, mean_idx].values.astype(np.float)

        # FIXME: Handle mis-aligned times here
        unique_names = np.unique(names)
        stack_times = None
        stack_areas = []
        stack_means = []
        for name in unique_names:
            mask = names == name
            if stack_times is None:
                stack_times = times[mask]
            elif not np.allclose(stack_times, times[mask]):
                raise ValueError(f'Misaligned time vectors in {infile}')
            stack_areas.append(areas[mask])
            stack_means.append(means[mask])

        # Now we have the times and area columns in the order we need them
        times = stack_times
        means = np.stack(stack_means, axis=1)
        areas = np.stack(stack_areas, axis=1)

        # Area == 0 means that the ROI isn't valid for this timepoint
        means[areas <= 0] = np.nan

        # Time < 0 means that this sample isn't valid for this timepoint
        means[times < 0] = np.nan

        keep_cols = np.sum(~np.isnan(means), axis=0) >= min_samples
        means = means[:, keep_cols]
        if means.shape[1] < 1:
            err = f'Empty stat file. No ROIs with at least {min_samples} samples in {infile}'
            raise OSError(err)
        return times, means

# Functions

def save_outfile(outfile: pathlib.Path,
                 time: np.ndarray,
                 signals: np.ndarray):
    """ Save the filtered output

    :param Path outfile:
        The path to the filtered output CSV or excel file
    :param ndarray time:
        The n x 1 array of time steps to save
    :param ndarray signals:
        The n x m array of signals to save
    """

    columns = ['Time']

    df = {'Time': time}
    for i in range(signals.shape[1]):
        key = f'Mean{i + 1:d}'
        columns.append(key)
        df[key] = signals[:, i]

    df = pd.DataFrame(df)
    if outfile.suffix in ('.csv', ):
        df.to_csv(str(outfile), columns=columns, index=False)
    elif outfile.suffix in ('.xls', '.xlsx'):
        df.to_excel(str(outfile), columns=columns, index=False)
    else:
        raise KeyError(f'Unknown output file suffix "{outfile}"')

def save_final_stats(outfile: pathlib.Path, stats: Dict):
    """ Save the composite stats for everything

    :param Path outfile:
        The path to the stat file to save
    :param dict stats:
        The dictionary of stats to save
    """

    df = {}
    for key, stat in sorted(stats.items()):
        if key == 'signal_stats':
            continue
        for i, rec in enumerate(stat):
            df.setdefault('Filename', []).append(key)
            df.setdefault('Trace', []).append(i + 1)
            for k, v in rec.items():
                df.setdefault(k, []).append(v)

    # Handle failed analyses
    if 'Filename' not in df:
        raise ValueError(f'No summary results from analyzing "{outfile}"')

    # Validate the final record numbers
    num_recs = len(df['Filename'])
    invalid_records = 0
    for key, vals in df.items():
        if len(vals) != num_recs:
            print(f'* Invalid {key} records: expected {num_recs}, got {len(vals)}')
            invalid_records += 1
    if invalid_records > 0:
        raise ValueError(f'Got {invalid_records} rows for "{outfile}"')

    # Save the final stats to a file
    print(f'Writing: {outfile}')
    df = pd.DataFrame(df)
    if outfile.suffix in ('.csv', ):
        df.to_csv(str(outfile), index=False)
    elif outfile.suffix in ('.xls', '.xlsx'):
        df.to_excel(str(outfile), index=False)
    else:
        raise KeyError(f'Unknown file type "{outfile}"')

def find_ca_data(datadir: pathlib.Path) -> Generator[pathlib.Path, None, None]:
    """ Find all the CA data in the data directory

    :param datadir:
        The data directory to search
    :returns:
        An iterator yielding paths to CA files
    """
    if datadir.is_file() and datadir.suffix in SUFFIXES:
        yield datadir
    else:
        for datafile in datadir.iterdir():
            if datafile.name.startswith('.'):
                continue
            if not datafile.is_file():
                continue
            # Ignore the outputs from the old analysis code
            if datafile.name == 'stats.xlsx':
                continue
            if datafile.suffix in SUFFIXES:
                yield datafile

max_stars_repo_name: david-a-joy/multilineage-organoid

max_issues_repo_name: david-a-joy/multilineage-organoid

max_forks_repo_name: david-a-joy/multilineage-organoid


text:
# -*- coding:utf-8 -*-
import os
import cv2
import time
import serial
import argparse
import numpy as np
from utils.anchor_generator import generate_anchors
from utils.anchor_decode import decode_bbox
from utils.nms import single_class_non_max_suppression
from load_model.pytorch_loader import load_pytorch_model, pytorch_inference

# model = load_pytorch_model('models/face_mask_detection.pth');
model = load_pytorch_model('./models/model360.pth')
# anchor configuration
feature_map_sizes = [[45, 45], [23, 23], [12, 12], [6, 6], [4, 4]]
anchor_sizes = [[0.04, 0.056], [0.08, 0.11], [0.16, 0.22], [0.32, 0.45], [0.64, 0.72]]
anchor_ratios = [[1, 0.62, 0.42]] * 5

# generate anchors
anchors = generate_anchors(feature_map_sizes, anchor_sizes, anchor_ratios)

# for inference , the batch size is 1, the model output shape is [1, N, 4],
# so we expand dim for anchors to [1, anchor_num, 4]
anchors_exp = np.expand_dims(anchors, axis=0)

id2class = {0: 'Mask', 1: 'NoMask'}

saved_img = False

def check_for_facemask(image,
              arduino,
              conf_thresh=0.5,
              iou_thresh=0.4,
              target_shape=(360, 360),
              draw_result=True):
    height, width, _ = image.shape
    image_resized = cv2.resize(image, target_shape)
    image_transposed = np.expand_dims(image_resized / 255.0, axis=0).transpose((0, 3, 1, 2))

    y_bboxes_output, y_cls_output = pytorch_inference(model, image_transposed)
    # remove the batch dimension, for batch is always 1 for inference.
    y_bboxes = decode_bbox(anchors_exp, y_bboxes_output)[0]
    y_cls = y_cls_output[0]
    # To speed up, do single class NMS, not multiple classes NMS.
    bbox_max_scores = np.max(y_cls, axis=1)
    bbox_max_score_classes = np.argmax(y_cls, axis=1)

    # keep_idx is the alive bounding box after nms.
    keep_idxs = single_class_non_max_suppression(y_bboxes,
                                                 bbox_max_scores,
                                                 conf_thresh=conf_thresh,
                                                 iou_thresh=iou_thresh)

    #if len(keep_idxs) == 0:
        #arduino.write(b'2')
        #time.sleep(1)
    for idx in keep_idxs:
        conf = float(bbox_max_scores[idx])
        class_id = bbox_max_score_classes[idx]
        bbox = y_bboxes[idx]
        # clip the coordinate, avoid the value exceed the image boundary.
        xmin = max(0, int(bbox[0] * width))
        ymin = max(0, int(bbox[1] * height))
        xmax = min(int(bbox[2] * width), width)
        ymax = min(int(bbox[3] * height), height)

        if draw_result:
            if class_id == 0:
                color = (0, 255, 0)
                #arduino.write(b'2')
            else:
                color = (255, 0, 0)
                #arduino.write(b'1')
                global saved_img
                if not saved_img:
                    dir = r'static\imgs'
                    cv2.imwrite(os.path.join(dir , str(len(os.listdir(dir))) + '.png' ), image)
                    saved_img = True
            #time.sleep(1)

            cv2.rectangle(image, (xmin, ymin), (xmax, ymax), color, 2)
            cv2.putText(image, "%s: %.2f" % (id2class[class_id], conf), (xmin + 2, ymin - 2),
                        cv2.FONT_HERSHEY_SIMPLEX, 0.8, color)

def run_on_video(video_path, conf_thresh):
    cap = cv2.VideoCapture(video_path)

    if not cap.isOpened():
        raise ValueError("Video open failed.")

    status = True
    while status:
        status, img_raw = cap.read()
        img_raw = cv2.cvtColor(img_raw, cv2.COLOR_BGR2RGB)
        if status:
            check_for_facemask(img_raw,
                      #arduino,
                      conf_thresh,
                      iou_thresh=0.5)
            cv2.imshow('image', img_raw[:, :, ::-1])
            cv2.waitKey(1)

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Face Mask Detection")
    parser.add_argument('--img-mode', type=int, default=1, help='set 1 to run on image, 0 to run on video.')
    parser.add_argument('--img-path', type=str, default='img/demo2.jpg', help='path to your image.')
    parser.add_argument('--video-path', type=str, default='0', help='path to your video, `0` means to use camera.')
    # parser.add_argument('--hdf5', type=str, help='keras hdf5 file')
    args = parser.parse_args()
    #arduino = serial.Serial('COM3', 9600, timeout = 0.1)
    #time.sleep(2)

    if args.img_mode:
        imgPath = args.img_path
        img = cv2.imread(imgPath)
        img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

        check_for_facemask(img, arduino)
    else:
        video_path = args.video_path
        if args.video_path == '0':
            video_path = 0
        run_on_video(video_path, conf_thresh=0.5)


text:
"""
Histogram Computation
=====================

This module provides a Job for the computation of the several histograms in
order to detect different changes in the scene (lightning, speaker motion, etc).

It also provides a Job for gathering all user input (slide location, speaker location).

.. autosummary::

  HistogramsLABDiff
  NumberOfVerticalStripesForSpeaker
  GenerateHistogramAreas

"""

from ..job import Job
import cv2
import numpy as np
import functools

from ....util.tools import get_polygon_outer_bounding_box, crop_image_from_normalized_coordinates, \
                           sort_dictionary_by_integer_key
from ....util.functor import Functor
from .select_polygon import SelectPolygonJob, SelectSlide, SelectSpeaker

class HistogramsLABDiff(Job):
    """
    Computes the histograms on the difference image of two consecutive frames.

    The difference image
    is expressed in the LAB color space using an Euclidean metric. The histograms are computed on
    several areas of the image plane.

    .. rubric:: Workflow inputs

    The inputs of the parents are:

    * A list of tuples `(name, rectangle)` specifying the locations where the histogram should
      be computed.

      * `name` is indicating the name of the rectangle.
      * `rectangle` is given as `(x,y, width, height)`, in **normalized coordinates**

      If several rectangles exist for the same name, those are merged (useful if the area
      is defined by several disconnected polygons). See :py:class:`GenerateHistogramAreas` for a possible
      input

    * A list of images specified by filename, on which the histograms will be computed

    .. rubric:: Workflow outputs

    The output of this Job is binary a function::

        frame_index, rectangle_name -> histogram

    that provides the histogram in the difference image in this particular rectangle.

    .. rubric:: Complexity

    Linear in the number of thumbnails. Reads each thumbnail image once.
    """

    #: Name of the job in the workflow
    name = 'histogram_imlabdiff'

    #: Cached outputs:
    #:
    #: * ``histograms_labdiff`` histogram of the difference image of two consecutive images in the video
    #:   sequence (computed in LAB space)
    outputs_to_cache = ['histograms_labdiff']

    def __init__(self,
                 *args,
                 **kwargs):
        super(HistogramsLABDiff, self).__init__(*args, **kwargs)

    def load_state(self):
        """
        Sort the histograms by ``frame_index`` in order to be able to compare states.

        This is necessary because the json module can load and store dictionaries
        out of order.
        """
        state = super(HistogramsLABDiff, self).load_state()

        if state is None:
            return None

        histograms_labdiff = state['histograms_labdiff']

        for area in histograms_labdiff.keys():
            histograms_labdiff[area] = sort_dictionary_by_integer_key(histograms_labdiff[area])

        state['histograms_labdiff'] = histograms_labdiff
        return state

    def run(self, *args, **kwargs):
        assert(len(args) >= 2)

        self.rectangle_locations = args[0]

        image_list = args[1]

        # init
        self.histograms_labdiff = {}

        rectangle_names = zip(*self.rectangle_locations)[0]
        unique_rectangle_names = list(set(rectangle_names))

        for name in unique_rectangle_names:
            element = self.histograms_labdiff.get(name, {})
            self.histograms_labdiff[name] = element

        # perform the computation
        im_index_tm1 = cv2.imread(image_list[0])
        imlab_index_tm1 = cv2.cvtColor(im_index_tm1, cv2.COLOR_BGR2LAB)

        for index, filename in enumerate(image_list[1:], 1):
            im_index_t = cv2.imread(filename)
            imlab_index_t = cv2.cvtColor(im_index_t, cv2.COLOR_BGR2LAB)

            # color diff
            im_diff = (imlab_index_t - imlab_index_tm1) ** 2
            im_diff_lab = np.sqrt(np.sum(im_diff, axis=2))

            # Compute histogram for every area
            for name, rect in self.rectangle_locations:
                cropped = crop_image_from_normalized_coordinates(im_diff_lab, rect)
                histogram = cv2.calcHist([cropped.astype(np.uint8)], [0], None, [256], [0, 256])

                # Merge histograms if necessary
                histogram_to_merge = self.histograms_labdiff[name].get(index, None)
                if histogram_to_merge is not None:
                    histogram += histogram_to_merge

                self.histograms_labdiff[name][index] = histogram

            # @note(Stephan):
            # The histograms are stored as a python list in order to serialize them via JSON.
            for name in unique_rectangle_names:
                histogram_np_array = self.histograms_labdiff[name][index]
                self.histograms_labdiff[name][index] = histogram_np_array.tolist()

    def get_outputs(self):
        super(HistogramsLABDiff, self).get_outputs()
        if self.histograms_labdiff is None:
            raise RuntimeError('The points have not been selected yet')

        return Functor(self.histograms_labdiff, transform=functools.partial(np.array, dtype=np.float32))

class NumberOfVerticalStripesForSpeaker(Job):
    """Indicates the number of vertical stripes used for speaker tracking.

    .. rubric:: Runtime parameters

    * ``nb_vertical_stripes`` the number of vertical stripes for the speaker tracking.

    """

    #: name of the job in the workflow
    name = 'number_of_vertical_stripes_speaker'

    #: Cached outputs:
    #:
    #: * ``nb_vertical_stripes`` number of vertical divisions of the image for speacker tracking
    outputs_to_cache = ['nb_vertical_stripes']

    def __init__(self, *args, **kwargs):
        super(NumberOfVerticalStripesForSpeaker, self).__init__(*args, **kwargs)
        assert('nb_vertical_stripes' in kwargs)

    def run(self, *args, **kwargs):
        pass

    def get_outputs(self):
        super(NumberOfVerticalStripesForSpeaker, self).get_outputs()

        return self.nb_vertical_stripes

class GenerateHistogramAreas(Job):
    """Generates the area of interest on which histograms should be computed by the downstream.

    .. rubric:: Workflow outputs

    Each area is identified by a name and a rectangle location.
    The output of this job is a list of tuples `(name, rect)` where each tuples
    contains:

    * The name of the area
    * A normalized rectangle `[x,y,width,height]` that specifies the area.
    """

    #: name of the job in the workflow
    name = 'gather_selections'

    parents = [SelectSlide, SelectSpeaker, NumberOfVerticalStripesForSpeaker]

    #: Cached outputs:
    #:
    #: * ``rectangle_locations`` name and location of the areas on which the histograms should be computed
    outputs_to_cache = ['rectangle_locations']

    def run(self, *args, **kwargs):
        self.rectangle_locations = []

        # slide location gives the position where to look for the illumination
        # changes detection
        slide_loc = args[0]
        slide_rec = get_polygon_outer_bounding_box(slide_loc)
        x, y, width, height = slide_rec

        # those two areas are the left and right side of the slide area
        first_light_change_area = [0, y, x, height]
        second_light_change_area = [x + width, y, 1 - (x + width), height]

        # Unicode names in order to compare to the json file
        self.rectangle_locations += [u'slides', first_light_change_area], \
                                    [u'slides', second_light_change_area]

        # speaker location is divided into vertical stripes on the full horizontal
        # extent
        speaker_loc = args[1]
        speaker_rec = get_polygon_outer_bounding_box(speaker_loc)
        _, y, _, height = speaker_rec
        nb_vertical_stripes = args[2]

        width_stripes = 1.0 / nb_vertical_stripes
        for i in range(nb_vertical_stripes - 1):
            x_start = width_stripes * i
            rect_stripe = [x_start, y, width_stripes, height]
            self.rectangle_locations += [u'speaker_%.2d' % i,
                                         rect_stripe],

        # final stripe adjusted a bit to avoid getting out the image plane
        rect_stripe = [1 - width_stripes, y, width_stripes, height]
        self.rectangle_locations += [u'speaker_%.2d' % (nb_vertical_stripes - 1),
                                     rect_stripe],

    def get_outputs(self):
        super(GenerateHistogramAreas, self).get_outputs()

        if self.rectangle_locations is None:
            raise RuntimeError('The Areas we want to compute the histograms on have not been computed yet.')

        return self.rectangle_locations

max_issues_repo_name: raffienficiaud/livius

max_forks_repo_name: raffienficiaud/livius


text:
"""Reciprocal distribution."""
import numpy

from ..baseclass import Dist
from ..operators.addition import Add

class Reciprocal(Dist):
    """
    Reciprocal distribution.

    Args:
        lower (float, Dist):
            Lower threshold of distribution. Must be smaller than ``upper``.
        upper (float, Dist):
            Upper threshold of distribution.

    Examples:
        >>> distribution = chaospy.Reciprocal(2, 4)
        >>> distribution
        Reciprocal(lower=2, upper=4)
        >>> q = numpy.linspace(0, 1, 5)
        >>> distribution.inv(q).round(4)
        array([2.    , 2.3784, 2.8284, 3.3636, 4.    ])
        >>> distribution.fwd(distribution.inv(q)).round(4)
        array([0.  , 0.25, 0.5 , 0.75, 1.  ])
        >>> distribution.pdf(distribution.inv(q)).round(4)
        array([0.7213, 0.6066, 0.5101, 0.4289, 0.3607])
        >>> distribution.sample(4).round(4)
        array([3.1462, 2.166 , 3.8645, 2.7937])
        >>> distribution.mom(1).round(4)
        7.8433
    """

    def __init__(self, lower=1, upper=2):
        self._repr = {"lower": lower, "upper": upper}
        Dist.__init__(self, lower=lower, upper=upper)

    def _pdf(self, x, lower, upper):
        return 1./(x*numpy.log(upper/lower))

    def _cdf(self, x, lower, upper):
        return numpy.log(x/lower)/numpy.log(upper/lower)

    def _ppf(self, q, lower, upper):
        return numpy.e**(q*numpy.log(upper/lower) + numpy.log(lower))

    def _lower(self, lower, upper):
        return lower

    def _upper(self, lower, upper):
        return upper

    def _mom(self, k, lower, upper):
        return ((upper*numpy.e**k-lower*numpy.e**k)/(numpy.log(upper/lower)*(k+(k == 0))))**(k != 0)

max_stars_repo_name: lblonk/chaospy


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

""" Module to handle World Ocean Atlas (WOA) climatology
"""

from os.path import expanduser
import re
from datetime import datetime

import numpy as np
from numpy import ma
import netCDF4
from scipy.interpolate import interp1d
# RectBivariateSpline
from scipy.interpolate import griddata

from .utils import dbsource
from .common import cropIndices

# ============================================================================
def woa_profile(var, d, lat, lon, depth, cfg):
    # Must improve here. This try make sense if fail because there isn't an
    #   etopo file, but if fail for another reason, like there is no lat,
    #   it will loose time trying from_dap.
    try:
        woa = woa_profile_from_file(var, d, lat, lon, depth, cfg)
    except:
        try:
            woa = woa_profile_from_dap(var, d, lat, lon, depth, cfg)
        except:
            print("Couldn't make woa_comparison of %s" % var)
            return

    return woa

def woa_profile_from_dap(var, d, lat, lon, depth, cfg):
    """
    Monthly Climatologic Mean and Standard Deviation from WOA,
    used either for temperature or salinity.

    INPUTS
        time: [day of the year]
        lat: [-90<lat<90]
        lon: [-180<lon<180]
        depth: [meters]

    Reads the WOA Monthly Climatology NetCDF file and
    returns the corresponding WOA values of salinity or temperature mean and
    standard deviation for the given time, lat, lon, depth.
    """
    if lon < 0:
        lon = lon+360

    url = cfg['url']

    doy = int(d.strftime('%j'))
    dataset = open_url(url)

    dn = (np.abs(doy-dataset['time'][:])).argmin()
    xn = (np.abs(lon-dataset['lon'][:])).argmin()
    yn = (np.abs(lat-dataset['lat'][:])).argmin()

    if re.match(r'temperature\d?$', var):
        mn = ma.masked_values(dataset.t_mn.t_mn[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.t_mn.attributes['_FillValue'])
        sd = ma.masked_values(dataset.t_sd.t_sd[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.t_sd.attributes['_FillValue'])
        # se = ma.masked_values(dataset.t_se.t_se[dn, :, yn, xn].reshape(
        #    dataset['depth'].shape[0]), dataset.t_se.attributes['_FillValue'])
        # Use this in the future. A minimum # of samples
        # dd = ma.masked_values(dataset.t_dd.t_dd[dn, :, yn, xn].reshape(
        #    dataset['depth'].shape[0]), dataset.t_dd.attributes['_FillValue'])
    elif re.match(r'salinity\d?$', var):
        mn = ma.masked_values(dataset.s_mn.s_mn[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.s_mn.attributes['_FillValue'])
        sd = ma.masked_values(dataset.s_sd.s_sd[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.s_sd.attributes['_FillValue'])
        # dd = ma.masked_values(dataset.s_dd.s_dd[dn, :, yn, xn].reshape(
        #    dataset['depth'].shape[0]), dataset.s_dd.attributes['_FillValue'])
    zwoa = ma.array(dataset.depth[:])

    ind = (depth <= zwoa.max()) & (depth >= zwoa.min())
    # Mean value profile
    f = interp1d(zwoa[~ma.getmaskarray(mn)].compressed(), mn.compressed())
    mn_interp = ma.masked_all(depth.shape)
    mn_interp[ind] = f(depth[ind])
    # The stdev profile
    f = interp1d(zwoa[~ma.getmaskarray(sd)].compressed(), sd.compressed())
    sd_interp = ma.masked_all(depth.shape)
    sd_interp[ind] = f(depth[ind])

    output = {'woa_an': mn_interp, 'woa_sd': sd_interp}

    return output

def woa_track_from_file(d, lat, lon, filename, varnames=None):
    """ Temporary solution: WOA for surface track
    """
    d = np.asanyarray(d)
    lat = np.asanyarray(lat)
    lon = np.asanyarray(lon)

    lon[lon < 0] += 360

    doy = np.array([int(dd.strftime('%j')) for dd in d])

    nc = netCDF4.Dataset(expanduser(filename), 'r')

    if varnames is None:
        varnames = {}
        for v in nc.variables.keys():
            if nc.variables[v].dimensions == \
                    (u'time', u'depth', u'lat', u'lon'):
                        varnames[v] = v

    output = {}
    for v in varnames:
        output[v] = []

    for d_n, lat_n, lon_n in zip(doy, lat, lon):
        # Get the nearest point. In the future interpolate.
        n_d = (np.abs(d_n - nc.variables['time'][:])).argmin()
        n_x = (np.abs(lon_n - nc.variables['lon'][:])).argmin()
        n_y = (np.abs(lat_n - nc.variables['lat'][:])).argmin()

        for v in varnames:
            output[v].append(nc.variables[varnames[v]][n_d, 0, n_y, n_x])

    for v in varnames:
        output[v] = ma.fix_invalid(output[v])

    return output

class WOA_URL(object):
    def __init__(self):
        pass

class WOA_var_nc(object):
    """
    Reads the WOA Monthly Climatology NetCDF file and
    returns the corresponding WOA values of salinity or temperature mean and
    standard deviation for the given time, lat, lon, depth.
    """
    def __init__(self, source):
        self.ncs = source

        self.load_dims(dims=['lat', 'lon', 'depth'])
        self.set_keys()

    def __getitem__(self, item):
        return self.data[item]

    def keys(self):
        return self.KEYS

    def load_dims(self, dims):
        self.dims = {}
        for d in dims:
            self.dims[d] = self.ncs[0][d][:]
            for nc in self.ncs[1:]:
                assert (self.dims[d] == nc[d][:]).all()

        self.dims['time'] = []
        mfrac = 365/12.
        for nc in self.ncs:
            assert nc.variables['time'].size == 1
            assert nc.variables['time'].units[:6] == 'months'
            self.dims['time'].append(mfrac * (nc['time'][0] % 12))
        self.dims['time'] = np.array(self.dims['time'])

    def set_keys(self):
        """
        """
        self.KEYS = []
        for v in self.ncs[0].variables.keys():
            if self.ncs[0].variables[v].dimensions == \
                    (u'time', u'depth', u'lat', u'lon'):
                        S = self.ncs[0].variables[v].shape
                        for nc in self.ncs[1:]:
                            assert v in nc.variables
                            assert nc.variables[v].shape == S
                        self.KEYS.append(v)

    def crop(self, doy, depth, lat, lon, var):
        """ Crop a subset of the dataset for each var

            Given doy, depth, lat and lon, it returns the smallest subset
              that still contains the requested coordinates inside it.

            It handels special cases like a region around greenwich and
            the international date line.

            Accepts 0 to 360 and -180 to 180 longitude reference.

            It extends time and longitude coordinates, so simplify the use
               of series. For example, a ship track can be requested with
               a longitude sequence like [352, 358, 364, 369, 380], and
               the equivalent for day of year above 365.
        """
        dims, idx = cropIndices(self.dims, lat, lon, depth, doy)
        subset = {}
        for v in var:
            subset[v] = ma.asanyarray([
                self.ncs[tnn][v][0, idx['zn'], idx['yn'], idx['xn']]
                for tnn in idx['tn']])
        return subset, dims

    def nearest(self, doy, depth, lat, lon, var):
        output = {}
        dims, idx = cropIndices(self.dims, lat, lon, depth, doy)
        for v in var:
            output[v] = ma.masked_all((doy.size, depth.size, lat.size,
                lon.size), dtype='f')
            for tn_out, t in enumerate(doy):
                tn_in = np.absolute(dims['time']-t).argmin()
                subset = self.ncs[tn_in][v][0, idx['zn'], idx['yn'], idx['xn']]
                for yn_out, y in enumerate(lat):
                    yn_in = np.absolute(dims['lat']-y).argmin()
                    for xn_out, x in enumerate(lon):
                        xn_in = np.absolute(dims['lon']-x).argmin()
                        for zn_out, z in enumerate(depth):
                            zn_in = np.absolute(dims['depth']-z).argmin()
                            output[v][tn_out, zn_out, yn_out, xn_out] = \
                                    subset[zn_in, yn_in, xn_in]
        return output

    def interpolate(self, doy, depth, lat, lon, var):
        """ Interpolate each var on the coordinates requested

        """
        subset, dims = self.crop(doy, depth, lat, lon, var)

        # Subset contains everything requested. No need to interpolate.
        if np.all([d in dims['time'] for d in doy]) & \
                np.all([z in dims['depth'] for z in depth]) & \
                np.all([y in dims['lat'] for y in lat]) & \
                np.all([x in dims['lon'] for x in lon]):
                    dn = np.nonzero([d in doy for d in dims['time']])[0]
                    zn = np.nonzero([z in depth for z in dims['depth']])[0]
                    yn = np.nonzero([y in lat for y in dims['lat']])[0]
                    xn = np.nonzero([x in lon for x in dims['lon']])[0]
                    output = {}
                    for v in subset:
                        # output[v] = subset[v][dn, zn, yn, xn]
                        # Seriously that this is the way to do it?!!??
                        output[v] = subset[v][:, :, :, xn][:, :, yn][:, zn][dn]
                    return output

        output = {}
        for v in var:
            output[v] = ma.masked_all(
                    (doy.size, depth.size, lat.size, lon.size),
                    dtype=subset[v].dtype)

            # These interpolators don't understand Masked Arrays, but do NaN
            if subset[v].dtype in ['int32']:
                subset[v] = subset[v].astype('f')
            subset[v][ma.getmaskarray(subset[v])] = np.nan
            subset[v] = subset[v].data

        # First linear interpolate on time.
        if not (doy == dims['time']).all():
            for v in subset.keys():
                f = interp1d(dims['time'], subset[v], axis=0)
                subset[v] = f(doy)
            dims['time'] = np.atleast_1d(doy)

        if not (np.all(lat == dims['lat']) and np.all(lon == dims['lon'])):
            # Lat x Lon target coordinates are the same for all time and depth.
            points_out = []
            for latn in lat:
                for lonn in lon:
                    points_out.append([latn, lonn])
            points_out = np.array(points_out)

            # Interpolate on X/Y plane
            for v in subset:
                tmp = np.nan * np.ones(
                        (doy.size, dims['depth'].size, lat.size, lon.size),
                        dtype=subset[v].dtype)
                for nt in range(doy.size):
                    for nz in range(dims['depth'].size):
                        data = subset[v][nt, nz]
                        # The valid data
                        idx = np.nonzero(~np.isnan(data))
                        if idx[0].size > 0:
                            points = np.array([
                                dims['lat'][idx[0]], dims['lon'][idx[1]]]).T
                            values = data[idx]
                            # Interpolate along the dimensions that have more than
                            #   one position, otherwise it means that the output
                            #   is exactly on that coordinate.
                            #ind = np.array([np.unique(points[:, i]).size > 1
                            #    for i in range(points.shape[1])])
                            #assert ind.any()

                            try:
                                values_out = griddata(
                                    #np.atleast_1d(np.squeeze(points[:, ind])),
                                    np.atleast_1d(np.squeeze(points)),
                                    values,
                                    #np.atleast_1d(np.squeeze(points_out[:, ind])))
                                    np.atleast_1d(np.squeeze(points_out)))
                            except:
                                values_out = []
                                for p in points_out:
                                    try:
                                        values_out.append(griddata(
                                            np.atleast_1d(np.squeeze(points)),
                                            values,
                                            np.atleast_1d(np.squeeze(
                                                p))))
                                    except:
                                        values_out.append(np.nan)
                                values_out = np.array(values_out)

                            # Remap the interpolated value back into a 4D array
                            idx = np.isfinite(values_out)
                            for [y, x], out in zip(
                                    points_out[idx], values_out[idx]):
                                tmp[nt, nz, y==lat, x==lon] = out
                subset[v] = tmp

        # Interpolate on z
        same_depth = (np.shape(depth) == dims['depth'].shape) and \
                np.allclose(depth, dims['depth'])
        if not same_depth:
            for v in list(subset.keys()):
                try:
                    f = interp1d(dims['depth'], subset[v], axis=1, bounds_error=False)
                    # interp1d does not handle Masked Arrays
                    subset[v] = f(np.array(depth))
                except:
                    print("Fail to interpolate '%s' in depth" % v)
                    del(subset[v])

        for v in subset:
            if output[v].dtype in ['int32']:
                subset[v] = np.round(subset[v])
            output[v][:] = ma.fix_invalid(subset[v][:])

        return output

    def extract(self, mode=None, **kwargs):
        """

            Possible scenarios:
              - Point:   doy{1},   depth{1},     lat{1},lon{1}
              - Profile: doy{1},   depth{0,1,n}, lat{1},lon{1}
              - Section: doy{1},   depth{0, n}, [lat{1},lon{n} | lat{n},lon{1}]

              - Track:   doy{1,n}, depth{1,n2},  lat{n},lon{n}
        """
        for k in kwargs:
            assert k in ['var', 'doy', 'depth', 'lat', 'lon'], \
                    "Wrong dimension to extract, check the manual"

        if 'var' in kwargs:
            var = np.atleast_1d(kwargs['var'])
        else:
            var = np.asanyarray(self.KEYS)

        doy = np.atleast_1d(kwargs['doy'])
        # This would only work if doy is 1D
        if type(doy[0]) is datetime:
            doy = np.array([int(d.strftime('%j')) for d in doy])

        if 'depth' in kwargs:
            depth = np.atleast_1d(kwargs['depth'])
        else:
            depth = self.dims['depth'][:]

        assert np.all(depth >= 0), "Depth was supposed to be positive."

        lat = np.atleast_1d(kwargs['lat'])
        lon = np.atleast_1d(kwargs['lon'])

        if mode == 'nearest':
            output = self.nearest(doy, depth, lat, lon, var)
        else:
            output = self.interpolate(doy, depth, lat, lon, var)

        for v in output:
                output[v] = np.atleast_1d(np.squeeze(output[v]))

        return output

    def track(self, mode=None, **kwargs):
        """

            Possible scenarios:
              - Track:   doy{1,n}, depth{1,n2}, lat{n}, lon{n}
        """
        for k in kwargs:
            assert k in ['var', 'doy', 'depth', 'lat', 'lon'], \
                    "Wrong dimension to extract, check the manual"

        if 'var' in kwargs:
            var = np.atleast_1d(kwargs['var'])
        else:
            var = np.asanyarray(self.KEYS)

        doy = np.atleast_1d(kwargs['doy'])
        if type(doy[0]) is datetime:
            doy = np.array([int(d.strftime('%j')) for d in doy])

        if 'depth' in kwargs:
            depth = np.atleast_1d(kwargs['depth'])
        else:
            depth = self.dims['depth'][:]

        assert np.all(depth >= 0), "Depth was supposed to be positive."

        lat = np.atleast_1d(kwargs['lat'])
        lon = np.atleast_1d(kwargs['lon'])

        assert lat.shape == lon.shape

        N = max(doy.size, depth.size, lat.size)

        if doy.shape == (1,):
            doy = doy * np.ones(N, dtype='i')
        if depth.shape == (1,):
            depth = depth * np.ones(N, dtype='i')
        if lat.shape == (1,):
            lat = lat * np.ones(N, dtype='i')
            lon = lon * np.ones(N, dtype='i')

        output = {}
        for v in var:
            output[v] = []

        for t, z, y, x in zip(doy, depth, lat, lon):
            if mode == 'nearest':
                tmp = self.nearest(
                        np.array([t]), np.array([z]), np.array([y]), np.array([x]), var)
            else:
                tmp = self.interpolate(
                        np.array([t]), np.array([z]), np.array([y]), np.array([x]), var)

            for v in tmp:
                output[v].append(tmp[v])

        for v in output:
            output[v] = np.atleast_1d(ma.array(output[v]).squeeze())

        return output

    def get_profile(var, doy, depth, lat, lon):
        print("get_profile is deprecated. You should migrate to extract()")
        return extract(var=var, doy=doy, depth=depth, lat=lat, lon=lon)

class WOA(object):
    """
    """
    def __init__(self, dbname='WOA18', resolution=None, tscale=None):
        self.dbname = dbname
        self.data = {'sea_water_temperature': None,
                'sea_water_salinity': None,
                'dissolved_oxygen': None,
                }
        self.resolution = resolution
        self.tscale = tscale

    def keys(self):
        return self.data.keys()

    def __getitem__(self, item):
        if item in ['TEMP', 'temperature']:
            return self['sea_water_temperature']
        elif item in ['PSAL', 'salinity']:
            return self['sea_water_salinity']
        elif item in ['DOXY']:
            return self['dissolved_oxygen']

        if self.data[item] is None:
            self.data[item] = WOA_var_nc(source=dbsource(
                self.dbname, item, self.resolution, self.tscale))
        return self.data[item]

max_stars_repo_name: petejan/oceansdb

max_issues_repo_name: petejan/oceansdb

max_forks_repo_name: petejan/oceansdb

reason: import numpy,from numpy,from scipy