text:
import tensorflow as tf
import numpy as np
from PIL import Image
import os

# On my machine, I have to use the CPU
# because occured some problem when i using the GPU,
# so you can comment the code below for better speed
os.environ["CUDA_VISIBLE_DEVICES"] = "-1"

img = np.invert(Image.open("1.png").convert('L')).ravel().reshape((1, 784))

def test_regression():
    tf.reset_default_graph()
    with tf.Session() as sess:
        saver = tf.train.import_meta_graph(
            os.path.join(os.path.dirname(__file__), "model/regression/regression-1000.meta"))
        saver.restore(sess, tf.train.latest_checkpoint('./model/regression/'))

        graph = tf.get_default_graph()

        # print all tensor name
        print([n.name for n in graph.as_graph_def().node])

        x = graph.get_tensor_by_name("regression/x:0")
        feed_dict = {x: img}

        regression = graph.get_tensor_by_name("y:0")

        print(sess.run(regression, feed_dict).flatten().tolist())

def test_cnn():
    tf.reset_default_graph()
    with tf.Session() as sess:
        saver = tf.train.import_meta_graph(os.path.join(os.path.dirname(__file__), "model/cnn/cnn-1000.meta"))
        saver.restore(sess, tf.train.latest_checkpoint('./model/cnn/'))

        graph = tf.get_default_graph()

        # print all tensor name
        print([n.name for n in graph.as_graph_def().node])

        x = graph.get_tensor_by_name("cnn/x:0")
        keep_prob = graph.get_tensor_by_name("keep_prob:0")
        feed_dict = {x: img, keep_prob: 1.0}

        cnn = graph.get_tensor_by_name("y_conv:0")
        print(sess.run(cnn, feed_dict).flatten().tolist())

if __name__ == '__main__':
    test_cnn()
    test_regression()


text:
import logging

import matplotlib
import numpy as np

from matplotlib.figure import Figure

from dbnd import task
from dbnd_test_scenarios.test_common.targets.target_test_base import TargetTestBase

logger = logging.getLogger(__name__)

matplotlib.use("Agg")

@task
def t_f_matplotlib_return():
    # type: ()-> Figure
    import matplotlib.pyplot as plt

    fig = plt.figure()
    ax1 = fig.add_subplot(2, 2, 1)
    ax1.hist(np.random.randn(100), bins=20, alpha=0.3)
    ax2 = fig.add_subplot(2, 2, 2)
    ax2.scatter(np.arange(30), np.arange(30) + 3 * np.random.randn(30))
    fig.add_subplot(2, 2, 3)
    return fig

class TestMatplotlibFigureOutputs(TargetTestBase):
    def test_matplotlib_return(self):
        t_f_matplotlib_return.dbnd_run()

ext: py

max_stars_repo_name: ipattarapong/dbnd

max_issues_repo_name: ipattarapong/dbnd

max_forks_repo_name: ipattarapong/dbnd

reason: import numpy


text:
'''
File: xy_field.py
Project: analyze
Created Date: 17/02/2021
Author: Shun Suzuki
-----
Last Modified: 24/02/2021
Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp)
-----
Copyright (c) 2021 Hapis Lab. All rights reserved.

'''

from shared import setup_pyplot, get_40kHz_amp, print_progress
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import mpl_toolkits.axes_grid1
import os
import glob
import re

def plot_acoustic_field_2d(axes, acoustic_pressures_2d, observe_x, observe_y, resolution, ticks_step, cmap='jet'):
    heatmap = axes.pcolor(acoustic_pressures_2d, cmap=cmap)
    x_label_num = int(math.floor((observe_x[1] - observe_x[0]) / ticks_step)) + 1
    y_label_num = int(math.floor((observe_y[1] - observe_y[0]) / ticks_step)) + 1
    x_labels = [observe_x[0] + ticks_step * i for i in range(x_label_num)]
    y_labels = [observe_y[0] + ticks_step * i for i in range(y_label_num)]
    x_ticks = [ticks_step / resolution * i for i in range(x_label_num)]
    y_ticks = [ticks_step / resolution * i for i in range(y_label_num)]
    axes.set_xticks(np.array(x_ticks) + 0.5, minor=False)
    axes.set_yticks(np.array(y_ticks) + 0.5, minor=False)
    axes.set_xticklabels(x_labels, minor=False)
    axes.set_yticklabels(y_labels, minor=False)

    return heatmap

def calc(data_path):
    cond = pd.read_csv(filepath_or_buffer=os.path.join(data_path, 'cond.txt'), sep=",", header=None)
    sample_rate = cond.at[0, 1]
    mV_per_Pa = cond.at[2, 1]
    dt = 1.0 / sample_rate

    p = re.compile(r'x([+-]?\d+\.?\d+?)y([+-]?\d+\.?\d+?)z([+-]?\d+\.?\d+?).csv')

    x_axis = []
    y_axis = []

    total = 0
    for filepath in glob.glob(os.path.join(data_path, '*')):
        m = p.match(filepath.split(os.path.sep)[-1])
        if m is None:
            continue

        x = float(m.group(1))
        y = float(m.group(2))

        x_axis.append(x)
        y_axis.append(y)
        total += 1

    x_axis = sorted(set(x_axis))
    y_axis = sorted(set(y_axis))

    rms = pd.DataFrame(index=y_axis, columns=x_axis)

    c = 0
    for filepath in glob.glob(os.path.join(data_path, '*')):
        m = p.match(filepath.split(os.path.sep)[-1])
        if m is None:
            continue

        x = float(m.group(1))
        y = float(m.group(2))

        df = pd.read_csv(filepath_or_buffer=filepath, sep=",")
        sound = df['  A Max [mV]']

        rms.at[y, x] = get_40kHz_amp(sound, dt) / mV_per_Pa / np.sqrt(2)
        c += 1
        print_progress(c, total)

    print()
    rms.to_csv('xy.csv')

def plot(plot_r):
    rms = pd.read_csv('xy.csv', index_col=0)
    resolution = float(rms.columns[1]) - float(rms.columns[0])

    nx = len(rms.columns)
    ny = len(rms.index)
    nx_min = int((nx - 1) / 2 - plot_r / resolution / 2)
    nx_max = int((nx - 1) / 2 + plot_r / resolution / 2) + 1
    ny_min = int((ny - 1) / 2 - plot_r / resolution / 2)
    ny_max = int((ny - 1) / 2 + plot_r / resolution / 2) + 1

    rms = rms.to_numpy().transpose().astype(np.float32)
    rms = rms[nx_min:nx_max, ny_min:ny_max]
    plot_xr = (-plot_r / 2, plot_r / 2)
    plot_yr = (-plot_r / 2, plot_r / 2)

    print('max [Pa]: ', rms.max())
    fig = plt.figure(figsize=(7, 6), dpi=DPI)
    axes = fig.add_subplot(111, aspect='equal')
    heat_map = plot_acoustic_field_2d(axes, rms, plot_xr, plot_yr, resolution, ticks_step=10.0)
    divider = mpl_toolkits.axes_grid1.make_axes_locatable(axes)
    plt.ylabel(r'$y\,[\mathrm{mm}]$', fontname='Arial', fontsize=18)
    plt.xlabel(r'$x\,[\mathrm{mm}]$', fontname='Arial', fontsize=18)
    cax = divider.append_axes('right', size='5%', pad='3%')
    cax.tick_params(labelsize=16)
    fig.colorbar(heat_map, cax=cax)
    cax.set_ylabel(r'$\mathrm{RMS\ of\ acoustic\ pressure}\,[\mathrm{Pa}]$', fontsize=18)

    plt.tight_layout()
    plt.savefig(os.path.join('plot', 'xy' + ext), bbox_inches='tight', pad_inches=0)

if __name__ == '__main__':
    os.makedirs('plot', exist_ok=True)
    setup_pyplot()

    DPI = 300
    ext = '.pdf'
    calc('./raw_data/xy')
    plot(80)


text:
import numpy as np
import tensorly
import copy

def CP_ALS(tensor, rank):
    #http://tensorly.org/stable/modules/generated/tensorly.decomposition.CP.html#tensorly.decomposition.CP
    return tensorly.decomposition.CP(tensor, rank)

def eig0(mat, it_time=100, tol=1e-15):
    """
    :param mat: 输入矩阵（实对称阵）
    :param it_time: 最大迭代步数
    :param tol: 收敛阈值
    :return lm: （绝对值）最大本征值
    :return v1: 最大本征向量
    """
    # 初始化向量
    v1 = np.random.randn(mat.shape[0],)
    v0 = copy.deepcopy(v1)
    lm = 1
    for n in range(it_time):  # 开始循环迭代
        v1 = mat.dot(v0)  # 计算v1 = M V0
        lm = np.linalg.norm(v1)  # 求本征值
        v1 /= lm  # 归一化v1
        # 判断收敛
        conv = np.linalg.norm(v1 - v0)
        if conv < tol:
            break
        else:
            v0 = copy.deepcopy(v1)
    return lm, v1

def svd0(mat, it_time=100, tol=1e-15):
    """
    Recursive algorithm to find the dominant singular value and vectors
    :param mat: input matrix (assume to be real)
    :param it_time: max iteration time
    :param tol: tolerance of error
    :return u: the dominant left singular vector
    :return s: the dominant singular value
    :return v: the dominant right singular vector
    """
    dim0, dim1 = mat.shape
    # 随机初始化奇异向量
    u, v = np.random.randn(dim0, ), np.random.randn(dim1, )
    # 归一化初始向量
    u, v = u/np.linalg.norm(u), v/np.linalg.norm(v)
    s = 1

    for t in range(it_time):
        # 更新v和s
        v1 = u.dot(mat)
        s1 = np.linalg.norm(v1)
        v1 /= s1
        # 更新u和s
        u1 = mat.dot(v1)
        s1 = np.linalg.norm(u1)
        u1 /= s1
        # 计算收敛程度
        conv = np.linalg.norm(u - u1) / dim0 + np.linalg.norm(v - v1) / dim1
        u, s, v = u1, s1, v1
        # 判断是否跳出循环
        if conv < tol:
            break
    return u, s, v

# rank-1 decomposition
def rank1decomp(x, it_time=100, tol=1e-15):
    """
    :param x: 待分解的张量
    :param it_time: 最大迭代步数
    :param tol: 迭代终止的阈值
    :return vs: 储存rank-1分解各个向量的list
    :return k: rank-1系数
    """
    ndim = x.ndim  # 读取张量x的阶数
    dims = x.shape  # 读取张量x各个指标的维数

    # 初始化vs中的各个向量并归一化
    vs = list()  # vs用以储存rank-1分解得到的各个向量
    for n in range(ndim):
        _v = np.random.randn(dims[n])
        vs.append(_v / np.linalg.norm(_v))
    k = 1

    for t in range(it_time):
        vs0 = copy.deepcopy(vs)  # 暂存各个向量以计算收敛情况
        for _ in range(ndim):
            # 收缩前(ndim-1)个向量，更新最后一个向量
            x1 = copy.deepcopy(x)
            for n in range(ndim-1):
                x1 = np.tensordot(x1, vs[n], [[0], [0]])
            # 归一化得到的向量，并更新常数k
            k = np.linalg.norm(x1)
            x1 /= k
            # 将最后一个向量放置到第0位置
            vs.pop()
            vs.insert(0, x1)
            # 将张量最后一个指标放置到第0位置
            x = x.transpose([ndim-1] + list(range(ndim-1)))
        # 计算收敛情况
        conv = np.linalg.norm(np.hstack(vs0) - np.hstack(vs))
        if conv < tol:
            break
    return vs, k

max_stars_repo_name: NieCha/TensorNetwork

max_issues_repo_name: NieCha/TensorNetwork

max_forks_repo_name: NieCha/TensorNetwork


text:
from itertools import chain, accumulate

import numpy as np

class DataMatrix():
    """A matrix of data, with names for features and samples, and weights.

    Arguments:
        data(numpy.ndarray):
            Initialization of `.data`, the actual data to store. Shape
            should be samples by features.
        samples(List[str]):
            Initialization of sample names (`.samples`). Should be same length
            as first dimension of data. Default is a list of empty strings.
        features(List[str]):
            Initialization of features names (`.features`). Should be same
            length as second dimension of data. Default is a list of empty
            strings.
        weights(ndarray):
            Initialization of `.weights`. Should be same shape as data. Default
            is a weight of one for every datapoint.

    Attributes:
        data(numpy.ndarray):
            Array with the actual data.
        weights(numpy.ndarray):
            Array with precision weights for the data, same shape as `data`.
        samples(List[str]):
            List with sample names.
        features(List[str]):
            List with feature names.
    """

    def __init__(self, data, samples=None, features=None, weights=None):
        if samples is None:
            samples = ["" for _ in range(data.shape[0])]
        if data.shape[0] != len(samples):
            raise Exception("Data dimension 0 should be the same as the "
                            "number of sample names")

        if features is None:
            features = ["" for _ in range(data.shape[1])]
        if data.shape[1] != len(features):
            raise Exception("Data dimension 1 should be the same as the "
                            "number of feature names")

        if weights is None:
            weights = np.ones_like(data)
        if weights.shape != data.shape:
            raise Exception("Dimensions of weights does not equal "
                            "dimensions of data")

        self.data = data
        self.samples = samples
        self.features = features
        self.weights = weights

    @property
    def dataW(self):
        d = self.data * self.weights
        d -= np.mean(d, 0, keepdims=True)
        return d

class StackedDataMatrix(DataMatrix):
    """A `.DataMatrix` containing multiple data types.

    Data and weights matrices are concatenated along the feature dimension (1,
    columns). Lists of feature names are concatenated. Sample names should be
    equal and in the same order for all data types. Data types of all data and
    weight matrices should match.

    Arguments:
        matrices(List[DataMatrix]):
            List of `.DataMatrix` to stack together.
        dt_names(List[str]):
            Names of data types, same order and length as matrices. Defaults
            to empty strings.

    Attributes:
        dt_names(List[str]):
            names of data types.
        dt_n_features(List[str]):
            number of features per data type.
        slices(List[slice]):
            slices to use for indexing `~.data`.
    """

    def __init__(self, matrices, dt_names=None):
        if dt_names is None:
            dt_names = ["" for _ in range(len(matrices))]

        samples = matrices[0].samples
        for m in matrices:
            for n1, n2 in zip(m.samples, samples):
                if n1 != n2:
                    raise Exception("Samples names don't match.")

        dtype = matrices[0].data.dtype
        for m in matrices:
            if m.data.dtype != dtype:
                raise Exception("dtypes of matrices don't match")

        dt_n_features = [len(m.features) for m in matrices]
        slices = self._compute_slices(dt_n_features)

        data = np.zeros((len(samples), sum(dt_n_features)), dtype=dtype)
        weights = np.ones((len(samples), sum(dt_n_features)), dtype=dtype)
        for m, sl in zip(matrices, slices):
            data[:, sl] = m.data
            weights[:, sl] = m.weights
        features = list(chain(*[m.features for m in matrices]))

        super().__init__(data, samples, features, weights)

        self.dt_names = dt_names
        self.n_dt = len(dt_names)
        self.slices = slices
        self.dt_n_features = dt_n_features

    def dt(self, idx, weighted=True) -> np.ndarray:
        """Index the data to get one datatype.

        Attributes:
            idx(int or str):
                Name or integer index of datatype.
        Returns:
            Numpy array with the data of one data type.
        """
        if isinstance(idx, str):
            idx = self.dt_names.index(idx)
        d = self.data[:, self.slices[idx]]
        if weighted:
            d = d * self.weights[:, self.slices[idx]]
            d -= np.mean(d, 0, keepdims=True)
        return d

    @staticmethod
    def _compute_slices(feature_lengths):
        return [slice(i_start, i_stop) for i_start, i_stop in
                zip(chain([0], accumulate(feature_lengths)),
                    accumulate(feature_lengths))]

    def to_netcdf(self, fn, feature_names=None, sample_dim='sample'):
        import netCDF4

        if feature_names is None:
            feature_names = [f'{dt}_feature' for dt in self.dt_names]

        with netCDF4.Dataset(fn, 'w') as ds:
            ds.createDimension(sample_dim, len(self.samples))
            samples_a = np.array(self.samples)
            samples = ds.createVariable(sample_dim, samples_a.dtype,
                                        sample_dim)
            for i, s in enumerate(self.samples):
                samples[i] = s

            ds.setncattr_string('data_types', self.dt_names)
            for dt_i, dt in enumerate(self.dt_names):
                ds.createDimension(feature_names[dt_i],
                                   self.dt_n_features[dt_i])
                fvar_a = np.array(self.features[self.slices[dt_i]])
                fvar = ds.createVariable(feature_names[dt_i], fvar_a.dtype,
                                         feature_names[dt_i])
                for i, v in enumerate(self.features[self.slices[dt_i]]):
                    fvar[i] = v

                var = ds.createVariable(dt, 'f8',
                                        (sample_dim, feature_names[dt_i]))
                var[:] = self.data[:, self.slices[dt_i]]

                wvar = ds.createVariable(f'{dt}_weights', 'f8',
                                         (sample_dim, feature_names[dt_i]))
                wvar[:] = self.weights[:, self.slices[dt_i]]

    @classmethod
    def from_netcdf(cls, fn, sample_dim='sample'):
        import netCDF4

        with netCDF4.Dataset(fn, 'r') as ds:
            dt_names = ds.getncattr('data_types')
            samples = list(ds[sample_dim])
            dms = dict()
            for dt in dt_names:
                assert(len(ds[dt].dimensions) == 2)
                assert(ds[dt].dimensions[0] == sample_dim)
                feature_dim = ds[dt].dimensions[1]
                dms[dt] = DataMatrix(
                    np.array(ds[dt]),
                    samples,
                    list(ds[feature_dim]),
                    np.array(ds[f'{dt}_weights']),
                )

        return cls([dms[dt] for dt in dt_names], dt_names)

max_stars_repo_name: NKI-CCB/funcsfa

max_issues_repo_name: NKI-CCB/funcsfa

max_forks_repo_name: NKI-CCB/funcsfa


text:
# Missions.py
#
# Created:  Feb 2016, E. Botero
# Modified:

# ----------------------------------------------------------------------
#   Imports
# ----------------------------------------------------------------------

import SUAVE
from SUAVE.Core import Units
import time

import numpy as np

# ----------------------------------------------------------------------
#   Define the Mission
# ----------------------------------------------------------------------

def setup(analyses,vehicle):

    # the mission container
    missions = SUAVE.Analyses.Mission.Mission.Container()
    missions.mission = mission(analyses,vehicle)

    return missions

def mission(analyses,vehicle):

    # ------------------------------------------------------------------
    #   Initialize the Mission
    # ------------------------------------------------------------------

    mission     = SUAVE.Analyses.Mission.Sequential_Segments()
    mission.tag = 'mission'

    mission.atmosphere  = SUAVE.Attributes.Atmospheres.Earth.US_Standard_1976()
    mission.planet      = SUAVE.Attributes.Planets.Earth()

    # unpack Segments module
    Segments = SUAVE.Analyses.Mission.Segments

    # base segment
    base_segment = Segments.Segment()
    base_segment.process.iterate.initials.initialize_battery = SUAVE.Methods.Missions.Segments.Common.Energy.initialize_battery

    # ------------------------------------------------------------------
    #   Cruise Segment: Constant Dynamic Pressure, Constant Altitude
    # ------------------------------------------------------------------

    segment = SUAVE.Analyses.Mission.Segments.Cruise.Constant_Dynamic_Pressure_Constant_Altitude(base_segment)
    segment.tag = "cruise"

    # connect vehicle configuration
    segment.analyses.extend(analyses.base)

    # segment attributes
    segment.state.numerics.number_control_points = 32
    segment.dynamic_pressure = 115.0 * Units.pascals
    segment.start_time       = time.strptime("Tue, Jun 21  11:00:00  2020", "%a, %b %d %H:%M:%S %Y",)
    segment.altitude         = 1000.0 * Units.feet
    segment.distance         = 1000.0 * Units.km
    segment.charge_ratio     = 1.0
    segment.latitude         = 37.4
    segment.longitude        = -122.15
    segment.battery_energy   = vehicle.base.networks.solar_low_fidelity.battery.max_energy
    segment.state.conditions.frames.wind.body_rotations[:,2] = 125.* Units.degrees

    mission.append_segment(segment)

    return mission


text:
import pandas as pd
import numpy as np
import requests
import joblib
from sklearn.neighbors import KDTree
from sklearn.preprocessing import StandardScaler

cultdf = pd.read_csv('data/cultdf.csv')

def get_culture(id):
    get = requests.get(f"https://api.citrics.io/jkekal6d6e5si3i2ld66d4dl/citydata/{id}").json()
    age_dist = list(get['Age Distribution'].values())
    p2010 = get['Population Growth']['2010']
    p2018 = get['Population Growth']['2018']
    p2017 = get['Population Growth']['2017']
    p2016 = get['Population Growth']['2016']
    if p2010 != 0 and p2018 != 0:
        popgrow = round(((p2018 - p2010) / p2010) * 100, 2)
    else:
        popgrow = 0
    diveristy = list(get['Ethnicity'].values())
    education = list(get['Educational Attainment'].values())
    birth = get['Recent Mothers']['Birth Rate']['Avg']

    res = age_dist + [popgrow, p2016, p2017 p2018] + diveristy + education + [birth]
    return(res)

scaler = StandardScaler()
scaler.fit(cultdf.values)
data = scaler.transform(cultdf.values)

tree = KDTree(data, leaf_size=10)
test = scaler.transform(np.array(get_culture(17927)).reshape(1,-1))

for i in tree.query(test, k=6)[1][0]:
    print(requests.get(f'https://api.citrics.io/jkekal6d6e5si3i2ld66d4dl/citydata/{i}').json()['name_with_com'], i)

joblib.dump(scaler, 'culture_scaler.joblib')
joblib.dump(tree, 'culture_model.joblib')

lscaler = joblib.load('culture_scaler.joblib')
ltree = joblib.load('culture_model.joblib')

test = lscaler.transform(np.array(get_culture(17927)).reshape(1,-1))

for i in ltree.query(test, k=6)[1][0]:
    print(requests.get(f'https://api.citrics.io/jkekal6d6e5si3i2ld66d4dl/citydata/{i}').json()['name_with_com'], i)

max_stars_repo_name: nephylum/city-data-comparison-ds

max_issues_repo_name: nephylum/city-data-comparison-ds

max_forks_repo_name: nephylum/city-data-comparison-ds


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Method Manager for Electrical Resistivity Tomography (ERT)"""

import os.path
import numpy as np
import matplotlib.pyplot as plt

import pygimli as pg
from pygimli.frameworks import MeshMethodManager
from .ertModelling import ERTModelling, ERTModellingReference
from .ert import createInversionMesh, createGeometricFactors, estimateError
from pygimli.utils import getSavePath

class ERTManager(MeshMethodManager):
    """ERT Manager.

    Method Manager for Electrical Resistivity Tomography (ERT)

    Todo
    ----
        * 3d
        * 3dtopo
        * complex on/off
        * closed geometry
        * transdim
        * singularity removal
        * ERT specific inversion options:
            * ...
    """

    def __init__(self, data=None, **kwargs):
        """Create ERT Manager instance.

        Parameters
        ----------
        data: :gimliapi:`GIMLI::DataContainerERT` | str
            You can initialize the Manager with data or give them a dataset
            when calling the inversion.

        Other Parameters
        ----------------
        * useBert: bool [True]
            Use Bert forward operator instead of the reference implementation.
        * sr: bool [True]
            Calculate with singularity removal technique.
            Recommended but needs the primary potential.
            For flat earth cases the primary potential will be calculated
            analytical. For domains with topography the primary potential
            will be calculated numerical using a p2 refined mesh or
            you provide primary potentials with setPrimPot.
        """
        self.useBert = kwargs.pop('useBert', True)
        self.sr = kwargs.pop('sr', True)

        super().__init__(data=data, **kwargs)
        self.inv.dataTrans = pg.trans.TransLogLU()

    def setSingularityRemoval(self, sr=True):
        """Turn singularity removal on or off."""
        self.reinitForwardOperator(sr=True)

    def createForwardOperator(self, **kwargs):
        """Create and choose forward operator."""
        verbose = kwargs.pop('verbose', False)
        self.useBert = kwargs.pop('useBert', self.useBert)
        self.sr = kwargs.pop('sr', self.sr)
        if self.useBert:
            pg.verbose('Create ERTModelling FOP')
            fop = ERTModelling(sr=self.sr, verbose=verbose)
        else:
            pg.verbose('Create ERTModellingReference FOP')
            fop = ERTModellingReference(**kwargs)

        return fop

    def load(self, fileName):
        """Load ERT data.

        Forwarded to :py:mod:`pygimli.physics.ert.load`

        Parameters
        ----------
        fileName: str
            Filename for the data.

        Returns
        -------
        data: :gimliapi:`GIMLI::DataContainerERT`
        """
        self.data = pg.physics.ert.load(fileName)
        return self.data

    def createMesh(self, data=None, **kwargs):
        """Create default inversion mesh.

        Forwarded to :py:mod:`pygimli.physics.ert.createInversionMesh`
        """
        d = data or self.data

        if d is None:
            pg.critical('Please provide a data file for mesh generation')

        return createInversionMesh(d, **kwargs)

    def setPrimPot(self, pot):
        """Set primary potential from external is not supported anymore."""
        pg.critical("Not implemented.")

    def simulate(self, mesh, scheme, res, **kwargs):
        """Simulate an ERT measurement.

        Perform the forward task for a given mesh, resistivity distribution &
        measuring scheme and return data (apparent resistivity) or potentials.

        For complex resistivity, the apparent resistivities is complex as well.

        The forward operator itself only calculates potential values for the
        electrodes in the given data scheme.
        To calculate apparent resistivities, geometric factors (k) are needed.
        If there are no values k in the DataContainerERT scheme, the function
        tries to calculate them, either analytically or numerically by using a
        p2-refined version of the given mesh.

        TODO
        ----
        * 2D + Complex + SR

        Args
        ----
        mesh : :gimliapi:`GIMLI::Mesh`
            2D or 3D Mesh to calculate for.

        res : float, array(mesh.cellCount()) | array(N, mesh.cellCount()) |
              list
            Resistivity distribution for the given mesh cells can be:
            . float for homogeneous resistivity (e.g. 1.0)
            . single array of length mesh.cellCount()
            . matrix of N resistivity distributions of length mesh.cellCount()
            . resistivity map as [[regionMarker0, res0],
                                  [regionMarker0, res1], ...]

        scheme : :gimliapi:`GIMLI::DataContainerERT`
            Data measurement scheme.

        Keyword Args
        ------------
        verbose: bool[False]
            Be verbose. Will override class settings.
        calcOnly: bool [False]
            Use fop.calculate instead of fop.response. Useful if you want
            to force the calculation of impedances for homogeneous models.
            No noise handling. Solution is put as token 'u' in the returned
            DataContainerERT.
        noiseLevel: float [0.0]
            add normally distributed noise based on
            scheme['err'] or on noiseLevel if error>0 is not contained
        noiseAbs: float [0.0]
            Absolute voltage error in V
        returnArray: bool [False]
            Returns an array of apparent resistivities instead of
            a DataContainerERT
        returnFields: bool [False]
            Returns a matrix of all potential values (per mesh nodes)
            for each injection electrodes.

        Returns
        -------
        DataContainerERT | array(data.size()) | array(N, data.size()) |
        array(N, mesh.nodeCount()):
            Data container with resulting apparent resistivity data and
            errors (if noiseLevel or noiseAbs is set).
            Optional returns a Matrix of rhoa values
            (for returnArray==True forces noiseLevel=0).
            In case of a complex valued resistivity model, phase values are
            returned in the DataContainerERT (see example below), or as an
            additionally returned array.

        Examples
        --------
        # >>> from pygimli.physics import ert
        # >>> import pygimli as pg
        # >>> import pygimli.meshtools as mt
        # >>> world = mt.createWorld(start=[-50, 0], end=[50, -50],
        # ...                        layers=[-1, -5], worldMarker=True)
        # >>> scheme = ert.createData(
        # ...     elecs=pg.utils.grange(start=-10, end=10, n=21),
        # ...     schemeName='dd')
        # >>> for pos in scheme.sensorPositions():
        # ...     _= world.createNode(pos)
        # ...     _= world.createNode(pos + [0.0, -0.1])
        # >>> mesh = mt.createMesh(world, quality=34)
        # >>> rhomap = [
        # ...     [1, 100. + 0j],
        # ...     [2, 50. + 0j],
        # ...     [3, 10.+ 0j],
        # ... ]
        # >>> data = ert.simulate(mesh, res=rhomap, scheme=scheme, verbose=1)
        # >>> rhoa = data.get('rhoa').array()
        # >>> phia = data.get('phia').array()
        """
        verbose = kwargs.pop('verbose', self.verbose)
        calcOnly = kwargs.pop('calcOnly', False)
        returnFields = kwargs.pop("returnFields", False)
        returnArray = kwargs.pop('returnArray', False)
        noiseLevel = kwargs.pop('noiseLevel', 0.0)
        noiseAbs = kwargs.pop('noiseAbs', 1e-4)
        seed = kwargs.pop('seed', None)
        sr = kwargs.pop('sr', self.sr)

        # segfaults with self.fop (test & fix)
        fop = self.createForwardOperator(useBert=self.useBert,
                                         sr=sr, verbose=verbose)
        fop.data = scheme
        fop.setMesh(mesh, ignoreRegionManager=True)

        rhoa = None
        phia = None

        isArrayData = False
        # parse the given res into mesh-cell-sized array
        if isinstance(res, int) or isinstance(res, float):
            res = np.ones(mesh.cellCount()) * float(res)
        elif isinstance(res, complex):
            res = np.ones(mesh.cellCount()) * res
        elif hasattr(res[0], '__iter__'):  # ndim == 2
            if len(res[0]) == 2:  # res seems to be a res map
                # check if there are markers in the mesh that are not defined
                # the rhomap. better signal here before it results in errors
                meshMarkers = list(set(mesh.cellMarkers()))
                mapMarkers = [m[0] for m in res]
                if any([mark not in mapMarkers for mark in meshMarkers]):
                    left = [m for m in meshMarkers if m not in mapMarkers]
                    pg.critical("Mesh contains markers without assigned "
                                "resistivities {}. Please fix given "
                                "rhomap.".format(left))
                res = pg.solver.parseArgToArray(res, mesh.cellCount(), mesh)
            else:  # probably nData x nCells array
                # better check for array data here
                isArrayData = True

        if isinstance(res[0], np.complex) or isinstance(res, pg.CVector):
            pg.info("Complex resistivity values found.")
            fop.setComplex(True)
        else:
            fop.setComplex(False)

        if not scheme.allNonZero('k') and not calcOnly:
            if verbose:
                pg.info('Calculate geometric factors.')
            scheme.set('k', fop.calcGeometricFactor(scheme))

        ret = pg.DataContainerERT(scheme)
        # just to be sure that we don't work with artifacts
        ret['u'] *= 0.0
        ret['i'] *= 0.0
        ret['r'] *= 0.0

        if isArrayData:
            rhoa = np.zeros((len(res), scheme.size()))
            for i, r in enumerate(res):
                rhoa[i] = fop.response(r)
                if verbose:
                    print(i, "/", len(res), " : ", pg.dur(), "s",
                          "min r:", min(r), "max r:", max(r),
                          "min r_a:", min(rhoa[i]), "max r_a:", max(rhoa[i]))
        else:  # res is single resistivity array
            if len(res) == mesh.cellCount():

                if calcOnly:
                    fop.mapERTModel(res, 0)

                    dMap = pg.core.DataMap()
                    fop.calculate(dMap)
                    if fop.complex():
                        pg.critical('Implement me')
                    else:
                        ret["u"] = dMap.data(scheme)
                        ret["i"] = np.ones(ret.size())

                    if returnFields:
                        return pg.Matrix(fop.solution())
                    return ret
                else:
                    if fop.complex():
                        res = pg.utils.squeezeComplex(res)

                    resp = fop.response(res)

                    if fop.complex():
                        rhoa, phia = pg.utils.toPolar(resp)
                    else:
                        rhoa = resp
            else:
                print(mesh)
                print("res: ", res)
                raise BaseException(
                    "Simulate called with wrong resistivity array.")

        if not isArrayData:
            ret['rhoa'] = rhoa

            if phia is not None:
                ret.set('phia', phia)
        else:
            ret.set('rhoa', rhoa[0])
            if phia is not None:
                ret.set('phia', phia[0])

        if returnFields:
            return pg.Matrix(fop.solution())

        if noiseLevel > 0:  # if errors in data noiseLevel=1 just triggers
            if not ret.allNonZero('err'):
                # 1A  and #100µV
                ret.set('err', self.estimateError(ret,
                                                  relativeError=noiseLevel,
                                                  absoluteUError=noiseAbs,
                                                  absoluteCurrent=1))
                print("Data error estimate (min:max) ",
                      min(ret('err')), ":", max(ret('err')))

            rhoa *= 1. + pg.randn(ret.size(), seed=seed) * ret('err')
            ret.set('rhoa', rhoa)

            ipError = None
            if phia is not None:
                if scheme.allNonZero('iperr'):
                    ipError = scheme('iperr')
                else:
                    # np.abs(self.data("phia") +TOLERANCE) * 1e-4absoluteError
                    if noiseLevel > 0.5:
                        noiseLevel /= 100.

                    if 'phiErr' in kwargs:
                        ipError = np.ones(ret.size()) * kwargs.pop('phiErr') \
                            / 1000
                    else:
                        ipError = abs(ret["phia"]) * noiseLevel

                    if verbose:
                        print("Data IP abs error estimate (min:max) ",
                              min(ipError), ":", max(ipError))

                phia += pg.randn(ret.size(), seed=seed) * ipError
                ret['iperr'] = ipError
                ret['phia'] = phia

        # check what needs to be setup and returned
        if returnArray:
            if phia is not None:
                return rhoa, phia
            else:
                return rhoa

        return ret

    def checkData(self, data=None):
        """Return data from container.

        THINKABOUT: Data will be changed, or should the manager keep a copy?
        """
        data = data or pg.DataContainerERT(self.data)
        if isinstance(data, pg.DataContainer):
            if not data.allNonZero('k'):
                pg.warn("Data file contains no geometric factors (token='k').")
                data['k'] = createGeometricFactors(data, verbose=True)

            if self.fop.complex():
                if not data.haveData('rhoa'):
                    pg.critical('Datacontainer have no "rhoa" values.')
                if not data.haveData('ip'):
                    pg.critical('Datacontainer have no "ip" values.')

                # pg.warn('check sign of phases')
                rhoa = data['rhoa']
                phia = -data['ip']/1000  # 'ip' is defined for neg mrad.
                # we should think about some 'phia' in rad

                return pg.utils.squeezeComplex(pg.utils.toComplex(rhoa, phia))

            else:
                if not data.haveData('rhoa'):
                    if data.allNonZero('r'):
                        pg.info("Creating apparent resistivies from "
                                "impedences rhoa = r * k")
                        data['rhoa'] = data['r'] * data['k']
                    elif data.allNonZero('u') and data.allNonZero('i'):
                        pg.info("Creating apparent resistivies from "
                                "voltage and currrent rhoa = u/i * k")
                        data['rhoa'] = data['u']/data['i'] * data['k']
                    else:
                        pg.critical("Datacontainer have neither: "
                                    "apparent resistivies 'rhoa', "
                                    "or impedances 'r', "
                                    "or voltage 'u' along with current 'i'.")

                if any(data['rhoa'] < 0) and \
                        isinstance(self.inv.dataTrans, pg.core.TransLog):
                    print(pg.find(data['rhoa'] < 0))
                    print(data['rhoa'][data['rhoa'] < 0])
                    pg.critical("Found negative apparent resistivities. "
                                "These can't be processed with logarithmic "
                                "data transformation. You should consider to "
                                "filter them out using "
                                "data.remove(data['rhoa'] < 0).")

                return data['rhoa']

        return data

    def checkErrors(self, err, dataVals):
        """Return relative error.

        Default we assume 'err' are relative vales.
        """
        if isinstance(err, pg.DataContainer):
            rae = None

            if not err.allNonZero('err'):
                pg.warn("Datacontainer have no 'err' values. "
                        "Fallback of 1mV + 3% using "
                        "ERTManager.estimateError(...) ")
                rae = self.estimateError(err, absoluteError=0.001,
                                         relativeError=0.03)
            else:
                rae = err['err']

            if self.fop.complex():

                ipe = None

                if err.haveData('iperr'):
                    amp, phi = pg.utils.toPolar(dataVals)
                    # assuming ipErr are absolute dPhi in mrad
                    ipe = err['iperr'] / abs((phi*1000))
                else:
                    pg.warn("Datacontainer have no 'iperr' values. "
                            "Fallback set to 0.01")
                    ipe = np.ones(err.size()) * 0.01

                return pg.cat(rae, ipe)

        return rae  # not set if err is no DataContainer (else missing)

    def estimateError(self, data=None, **kwargs):
        """Estimate error composed of an absolute and a relative part.

        Parameters
        ----------
        absoluteError : float [0.001]
            Absolute data error in Ohm m. Need 'rhoa' values in data.

        relativeError : float [0.03]
            relative error level in %/100

        absoluteUError : float [0.001]
            Absolute potential error in V. Need 'u' values in data. Or
            calculate them from 'rhoa', 'k' and absoluteCurrent if no 'i'
            is given

        absoluteCurrent : float [0.1]
            Current level in A for reconstruction for absolute potential V

        Returns
        -------
        error : Array
        """
        if data is None:  #
            error = estimateError(self.data, **kwargs)
            self.data["err"] = error
        else:  # the old way: better use ert.estimateError directly
            error = estimateError(data, **kwargs)

        return error

    def coverage(self):
        """Coverage vector considering the logarithmic transformation."""
        covTrans = pg.core.coverageDCtrans(self.fop.jacobian(),
                                           1.0 / self.inv.response,
                                           1.0 / self.inv.model)

        paramSizes = np.zeros(len(self.inv.model))
        for c in self.fop.paraDomain.cells():
            paramSizes[c.marker()] += c.size()

        return np.log10(covTrans / paramSizes)

    def standardizedCoverage(self, threshhold=0.01):
        """Return standardized coverage vector (0|1) using thresholding."""
        return 1.0*(abs(self.coverage()) > threshhold)

    def saveResult(self, folder=None, size=(16, 10), **kwargs):
        """Save all results in the specified folder.

        Saved items are:
            Inverted profile
            Resistivity vector
            Coverage vector
            Standardized coverage vector
            Mesh (bms and vtk with results)
        """
        subfolder = self.__class__.__name__
        path = getSavePath(folder, subfolder)

        pg.info('Saving resistivity data to: {}'.format(path))

        np.savetxt(path + '/resistivity.vector',
                   self.model)
        np.savetxt(path + '/resistivity-cov.vector',
                   self.coverage())
        np.savetxt(path + '/resistivity-scov.vector',
                   self.standardizedCoverage())

        m = pg.Mesh(self.paraDomain)
        m['Resistivity'] = self.paraModel(self.model)
        m['Resistivity (log10)'] = np.log10(m['Resistivity'])
        m['Coverage'] = self.coverage()
        m['S_Coverage'] = self.standardizedCoverage()
        m.exportVTK(os.path.join(path, 'resistivity'))
        m.saveBinaryV2(os.path.join(path, 'resistivity-pd'))
        self.fop.mesh().save(os.path.join(path, 'resistivity-mesh'))

        if self.paraDomain.dim() == 2:
            fig, ax = plt.subplots(figsize=size)
            self.showResult(ax=ax, coverage=self.coverage(), **kwargs)
            fig.savefig(path + '/resistivity.pdf', bbox_inches="tight")
            return path, fig, ax
        return path

if __name__ == "__main__":
    pass

max_stars_repo_name: JuliusHen/gimli

max_issues_repo_name: JuliusHen/gimli

max_forks_repo_name: JuliusHen/gimli


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Nov 27 19:46:57 2020

Grabs the cell parameters (orthogonal cells) from a LAMMPS trajectory.

Command-line options required:
    - file_in  : this is the trajectory
    - file_out : this is where the cell gets written

@author: andreypoletaev
"""

from sys import argv
from numpy import array, eye, savetxt

# =============================================================================
# %% read the first few lines of a trajectory and save the cell dimensions
# =============================================================================

## Parse inputs. Format: key=value
options = dict([ (x.split('=')[0], x.split('=')[1]) for x in argv[1:] ])
keys = list(options.keys())

assert 'file_in'  in keys, 'pass file_in=...  [path] as a command-line option'
assert 'file_out' in keys, 'pass file_out=... [path] as a command-line option'

with open(options['file_in'], 'r') as fin:

    for i in range(5) : l = fin.readline()

    cell = eye(3)

    for i in range(3) :
        l = fin.readline()
        d = array(l[:-1].split(' ')).astype(float)
        cell[i,i] = max(d) - min(d)

    savetxt(options['file_out'], cell, delimiter=' ', fmt='%.7g')

reason: from numpy


text:
import numpy as np
import scipy as sp

class ZCA():

    def __init__(self, bias=0.1):
        self.bias = bias

    def _flat(self, x):
        return np.reshape(x, (x.shape[0], -1))

    def _normalize(self, x):
        return (x.astype(np.float_) - self._mean) / self._std

    def fit(self, x):
        x = self._flat(x)
        self._mean = np.mean(x, axis=0, dtype=np.float64).astype(x.dtype)
        self._std = np.std(x, axis=0, dtype=np.float64).astype(x.dtype)
        x = self._normalize(x)
        try:
            # Perform dot product on GPU
            import cudarray as ca
            x_ca = ca.array(x)
            cov = np.array(ca.dot(x_ca.T, x_ca)).astype(np.float_)
        except:
            cov = np.dot(x.T, x)
        cov = cov / x.shape[0] + self.bias * np.identity(x.shape[1])
        s, v = np.linalg.eigh(cov)
        s = np.diag(1.0 / np.sqrt(s))
        self.whitener = np.dot(np.dot(v, s), v.T)
        return self

    def transform(self, x):
        shape = x.shape
        x = self._flat(x)
        x = self._normalize(x)
        x_white = np.dot(x, self.whitener.T)
        return np.reshape(x_white, shape)

class LCN():

    @staticmethod
    def gaussian_kernel(sigma, size=None):
        if size is None:
            size = int(np.ceil(sigma * 2.))
            if size % 2 == 0:
                size += 1
        xs = np.linspace(-size / 2., size / 2., size)
        kernel = 1 / (np.sqrt(2 * np.pi)) * \
            np.exp(-xs**2 / (2 * sigma**2)) / sigma
        return kernel / np.sum(kernel)

    def __init__(self, sigma, eps=1e-1, subtractive=False):
        self.sigma = sigma
        self.eps = eps
        self.subtractive = subtractive
        self.kernel = LCN.gaussian_kernel(sigma, size=int(sigma * 3))
        self.kernel = np.outer(self.kernel, self.kernel)

    def _normalize(self, x):
        return (x.astype(np.float_) - self._mean) / self._std

    def fit(self, x):
        self._mean = np.mean(x, dtype=np.float64).astype(x.dtype)
        self._std = np.std(x, dtype=np.float64).astype(x.dtype)
        return self

    def _transform_img(self, img):
        if img.ndim == 2:
            img = img[np.newaxis, :, :]
        n_channels = img.shape[0]

        # Calculate local mean
        mean = np.zeros((1,) + img.shape[1:])
        for i in range(n_channels):
            mean += sp.ndimage.filters.convolve(
                img[i], self.kernel, mode='nearest')
        mean /= n_channels

        # Center input with local mean
        centered = img - mean
        if self.subtractive:
            return centered

        # Calculate local standard deviation
        centered_sqr = centered**2
        std = np.zeros((1,) + img.shape[1:])
        for i in range(n_channels):
            std += sp.ndimage.filters.convolve(
                centered_sqr[i], self.kernel, mode='nearest')
        std /= n_channels
        std = np.sqrt(std)

        # Scale centered input with standard deviation
        return centered / (std + self.eps)

    def transform(self, x):
        x = self._normalize(x)
        return np.array([self._transform_img(img) for img in x])

reason: import numpy,import scipy


text:
# Copyright (c) 2020 - for information on the respective copyright owner
# see the NOTICE file and/or the repository https://github.com/boschresearch/blackboxopt
#
# SPDX-License-Identifier: Apache-2.0

from typing import Tuple

import numpy as np
import pytest

from blackboxopt import Evaluation, Objective
from blackboxopt.optimizers.staged.configuration_sampler import (
    StagedIterationConfigurationSampler,
)
from blackboxopt.optimizers.staged.iteration import StagedIteration
from blackboxopt.optimizers.staged.utils import greedy_promotion

class DummyOptimizer:
    def __init__(self, expected_configs_in_promote):
        self.expected_configs_in_promote = expected_configs_in_promote

    def decide_promotion(self, data, num_configs):
        assert len(data) == self.expected_configs_in_promote
        losses = [d.loss for d in data]
        ranks = np.argsort(np.argsort(losses))

        n = min(num_configs, len(data))
        return [datum.config_key for rank, datum in zip(ranks, data) if rank < n]

class DummyConfigurationSampler(StagedIterationConfigurationSampler):
    def __init__(self):
        self.num_configs = 0

    def sample_configuration(self) -> Tuple[dict, dict]:
        self.num_configs += 1
        return {"value": self.num_configs - 1}, {}

    def digest_evaluation(self, evaluation: Evaluation):
        """This dummy doesn't do anything with evaluation results."""
        pass

@pytest.mark.timeout(1)
def test_staged_iteration_get_and_digest_configuration():

    n_eval_specs = 3
    iteration = StagedIteration(
        iteration=0,
        num_configs=[n_eval_specs, 1],
        fidelities=[0.1, 1],
        config_sampler=DummyConfigurationSampler(),
        config_promotion_function=greedy_promotion,
        objective=Objective("loss", False),
    )

    eval_specs = []
    for i in range(n_eval_specs - 1):
        eval_specs.append(iteration.generate_evaluation_specification())
        assert eval_specs[-1].configuration["value"] == i
        iteration.digest_evaluation(
            eval_specs[-1].optimizer_info["id"],
            eval_specs[-1].create_evaluation(objectives={"loss": i}),
        )
        # make sure that the stage has not completed yet
        assert iteration.current_stage == 0

    # get the last config in this stage
    i = n_eval_specs - 1
    eval_specs.append(iteration.generate_evaluation_specification())
    assert eval_specs[-1].configuration["value"] == i
    assert iteration.current_stage == 0
    # make sure no configuration is returned when all configurations in a stage have
    # been queried, but not all have finished
    assert iteration.generate_evaluation_specification() is None

    # digest it and see that the next stage is reached
    iteration.digest_evaluation(
        eval_specs[-1].optimizer_info["id"],
        eval_specs[-1].create_evaluation(objectives={"loss": i}),
    )
    assert iteration.current_stage == 1
    assert not iteration.finished

    final_eval_spec = iteration.generate_evaluation_specification()
    assert final_eval_spec.configuration["value"] == 0
    iteration.digest_evaluation(
        final_eval_spec.optimizer_info["id"],
        final_eval_spec.create_evaluation(objectives={"loss": 0.0}),
    )
    assert iteration.finished
    assert iteration.generate_evaluation_specification() is None

@pytest.mark.timeout(1)
def test_staged_iteration_get_and_digest_configuration_with_crashes():

    n_eval_specs = 3
    opt = DummyOptimizer(n_eval_specs - 1)
    iteration = StagedIteration(
        iteration=0,
        num_configs=[n_eval_specs, n_eval_specs],
        fidelities=[1, 1],
        config_sampler=DummyConfigurationSampler(),
        config_promotion_function=opt.decide_promotion,
        objective=Objective("loss", False),
    )

    eval_specs = []
    for i in range(n_eval_specs):
        eval_specs.append(iteration.generate_evaluation_specification())
        assert eval_specs[-1].configuration["value"] == i

    for i in range(n_eval_specs - 1):
        iteration.digest_evaluation(
            eval_specs[i].optimizer_info["id"],
            eval_specs[i].create_evaluation(objectives={"loss": i}),
        )

    iteration.digest_evaluation(
        eval_specs[n_eval_specs - 1].optimizer_info["id"],
        eval_specs[n_eval_specs - 1].create_evaluation(objectives={"loss": None}),
    )

    assert all([e.status == "FINISHED" for e in iteration.evaluation_data[0][:-1]])
    assert iteration.evaluation_data[0][-1].status == "CRASHED"
    assert iteration.current_stage == 1

    # make sure the crashed config is not promoted and a new one is sampled instead
    eval_specs = []
    for i in range(n_eval_specs):
        eval_specs.append(iteration.generate_evaluation_specification())

    assert all([c.configuration["value"] == i for i, c in enumerate(eval_specs[:-1])])
    assert eval_specs[-1].configuration["value"] == n_eval_specs

max_stars_repo_name: boschresearch/blackboxopt

max_issues_repo_name: boschresearch/blackboxopt

max_forks_repo_name: boschresearch/blackboxopt


text:
import numpy as np
import pandas as pd
from sklearn.datasets import make_blobs
from gam import gam

def test_find_optimal_clusters():
    """"Create sample attributions with 4 clusters"""
    attributions_file = 'test_opt.csv'
    X, y = make_blobs(n_samples=40, n_features=2, centers=4,
                      cluster_std=0.01, center_box=(0.0, 1.0),
                      shuffle=True, random_state=42)

    print('blobs made - ', X.shape, y.shape)
    df = pd.DataFrame(columns=['x1', 'x2'], data=X)
    df.to_csv(attributions_file)

    """"Check cluster search via silhouette score"""
    g = gam.GAM(attributions_path=attributions_file, distance="kendall_tau")
    g.get_optimal_clustering()

    print('attributions file - ', g.attributions_path)
    print('data size = ', g.normalized_attributions.shape)
    print('what we settled on - ', g.k)
    assert(g.k == 4)

max_stars_repo_name: ceenam/global-attribution-mapping

max_issues_repo_name: ceenam/global-attribution-mapping

max_forks_repo_name: ceenam/global-attribution-mapping


text:
from __future__ import absolute_import, print_function

import numpy as np

import os

from PyQt5 import Qt, QtCore, QtGui, QtWidgets

import logging

logger = logging.getLogger(__name__)

class FloatSlider(QtWidgets.QSlider):
    floatValueChanged = QtCore.pyqtSignal(float)
    def __init__(self,*args):
        super(FloatSlider,self).__init__(*args)
        self.setRange(0,1.,100)
        self.valueChanged.connect(self.onChanged)

        self.setStyleSheet("QToolTip { color:white;}")

    def setRange(self,minVal,maxVal,steps=100):
        assert(minVal <= maxVal)
        super(FloatSlider,self).setRange(0,steps)
        self.minVal = minVal
        self.maxVal = maxVal
        self.steps = steps

    def _from_float(self,x):
        ind = int(self.steps*(x-self.minVal)/(self.maxVal-self.minVal))
        ind = max(0,min(self.steps,ind))
        logger.debug("floatslider (id = %s):  index from float %s = %s"%(id(self), x,ind))
        return ind

    def _from_int(self,n):
        return self.minVal+1.*(self.maxVal-self.minVal)*n/self.steps

    def setValue(self,val):
        logger.debug("floatslider (id = %s): setValue to : %s"%(id(self),val))
        self.floatValue = val

        # this surpressed anoying feedback
        old = self.blockSignals(True)
        super(FloatSlider,self).setValue(self._from_float(val))
        self.blockSignals(old)

    def value(self):
        return self.floatValue
        # return self._from_int(super(FloatSlider,self).value())

    def onChanged(self,ind):
        self.floatValue = self._from_int(ind)
        self.floatValueChanged.emit(self.floatValue)
        # self.floatValueChanged.emit(self._from_int(ind))

max_stars_repo_name: funkelab/spimagine

max_issues_repo_name: funkelab/spimagine

max_forks_repo_name: funkelab/spimagine


text:
import pandas as pd
from astropy.io import fits
import numpy as np
import tqdm.auto as tqdm
from collections import OrderedDict
from pathlib import Path

def dict_from_header(filename):
    """
    Parse a FITS header and extract the keys and values as an ordered dictionary. Multi-line keys like ``COMMENTS`` and ``HISTORY`` will be combined with commas. The resolved path will be inserted with the ``path`` key.

    Parameters
    ----------
    filename : str
        FITS file to parse

    Returns
    -------
    OrderedDict
    """
    summary = OrderedDict()
    summary["path"] = Path(filename).resolve()

    header = fits.getheader(filename)
    multi_entry_keys = {"COMMENT": [], "HISTORY": []}
    for k, v in header.items():
        if k == "":
            continue
        if k in multi_entry_keys:
            multi_entry_keys[k].append(v.lstrip())
        summary[k] = v

    for k, l in multi_entry_keys.items():
        if len(l) > 0:
            summary[k] = ", ".join(l)

    return summary

def observation_table(filenames, **kwargs):
    rows = [dict_from_header(filename) for filename in filenames]
    return pd.DataFrame(rows)

def parallactic_angle(header):
    if "D_IMRPAD" in header:
        return header["D_IMRPAD"] + header["LONPOLE"] - header["D_IMRPAP"]
    else:
        return parallactic_angle_altaz(header["ALTITUDE"], header["AZIMUTH"])

def parallactic_angle_hadec(ha, dec, lat=19.823806):
    """
    Calculate parallactic angle using the hour-angle and declination directly

    .. math::

        \\theta_\\mathrm{PA} = \\atan2{\\frac{\\sin\\theta_\\mathrm{HA}}{\\tan\\theta_\mathrm{lat}\\cos\\delta - \\sin\\delta \\cos\\theta_\\mathrm{HA}}}

    Parameters
    ----------
    ha : float
        hour-angle, in hour angles
    dec : float
        declination in degrees
    lat : float, optional
        latitude of observation in degrees, by default 19.823806

    Returns
    -------
    float
        parallactic angle, in degrees East of North
    """
    _ha = ha * np.pi / 12  # hour angle to radian
    _dec = np.deg2rad(dec)
    _lat = np.deg2rad(lat)
    sin_ha, cos_ha = np.sin(_ha), np.cos(_ha)
    sin_dec, cos_dec = np.sin(_dec), np.cos(_dec)
    pa = np.arctan2(sin_ha, np.tan(_lat) * cos_dec - sin_dec * cos_ha)
    return np.rad2deg(pa)

def parallactic_angle_altaz(alt, az, lat=19.823806):
    """
    Calculate parallactic angle using the altitude/elevation and azimuth directly

    .. math::
        \\theta_\\mathrm{PA} = \\atan2{\\frac{\\sin\\theta_\\mathrm{HA}}{\\tan\\theta_\mathrm{lat}\\cos\\delta - \\sin\\delta \\cos\\theta_\\mathrm{HA}}}

    Parameters
    ----------
    alt : float
        altitude or elevation, in degrees
    az : float
        azimuth, in degrees CCW from North
    lat : float, optional
        latitude of observation in degrees, by default 19.823806

    Returns
    -------
    float
        parallactic angle, in degrees East of North
    """
    ## Astronomical Algorithms, Jean Meeus
    # get angles, rotate az to S
    _az = np.deg2rad(az) - np.pi
    _alt = np.deg2rad(alt)
    _lat = np.deg2rad(lat)
    # calculate values ahead of time
    sin_az, cos_az = np.sin(_az), np.cos(_az)
    sin_alt, cos_alt = np.sin(_alt), np.cos(_alt)
    sin_lat, cos_lat = np.sin(_lat), np.cos(_lat)
    # get declination
    dec = np.arcsin(sin_alt * sin_lat - cos_alt * cos_lat * cos_az)
    # get hour angle
    ha = np.arctan2(sin_az, cos_az * sin_lat + np.tan(_alt) * cos_lat)
    # get parallactic angle
    pa = np.arctan2(np.sin(ha), np.tan(_lat) * np.cos(dec) - np.sin(dec) * np.cos(ha))
    return np.rad2deg(pa)

max_stars_repo_name: mileslucas/vampires-dpp

max_issues_repo_name: mileslucas/vampires-dpp

max_forks_repo_name: mileslucas/vampires-dpp

reason: import numpy,from astropy


text:
from task2vec import Task2Vec
from models import get_model
import datasets
import task_similarity
import argparse
import torch
import torch.nn.functional as F
import torchvision.transforms as tt
import torch.nn as nn
from tqdm import tqdm
from synbols_utils import Synbols
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import datetime

class HardAttentionRaw(nn.Module):
    """Attention through Raw Params
            Input: Param shape, No of Params"""
    def __init__(self, shape):
        super(AttentionRaw, self).__init__()
        # Attention Params: Shape N_attn, M as N Attentions of dimension M
        self.params = nn.Parameter(torch.randn(N, task2vec_dim), requires_grad = True)

    def forward(self, x):
        # Weighted/Atttention on x using Attention 'i'
        # X shape: N_z, M, Params shape: N_attn, M
        # Do Element wise multiplication
        return torch.mul(torch.unsqueeze(self.params, 1), x)

class AttentionRaw(nn.Module):
    """Attention through Raw Params
            Input: Param shape, No of Params"""
    def __init__(self, shape):
        super(AttentionRaw, self).__init__()
        # Attention Params: Shape N_attn, M as N Attentions of dimension M
        self.params = nn.Parameter(torch.randn(N, task2vec_dim), requires_grad = True)

    def forward(self, x):
        # Weighted/Atttention on x using Attention 'i'
        # X shape: N_z, M, Params shape: N_attn, M
        # Do Element wise multiplication
        soft = F.softmax(self.params, dim = -1)
        return torch.mul(torch.unsqueeze(soft, 1), x)

def cos_sim(A,B, eps=1e-8, ):
    A_norm = A.norm(dim = -1, keepdim = True)
    B_norm = B.norm(dim = -1, keepdim = True)
    A_norm = torch.maximum(A_norm, eps * torch.ones(A_norm.shape))
    B_norm = torch.maximum(B_norm, eps * torch.ones(B_norm.shape))
    A = torch.div(A,A_norm)
    B = torch.div(B,B_norm)
    return torch.mm(A,B.T)
def norm_unitvectors(M, eps = 1e-8):
    M_norm = M.norm(dim = -1, keepdim = True)
    M_norm = torch.maximum(M_norm, eps * torch.ones(M_norm.shape))
    return torch.div(M, M_norm)
def positives(Attrs):
    M = norm_unitvectors(Attrs)
    M_trans = M.transpose(1,2)
    return torch.bmm(M, M_trans)
def positive_loss(Attrs, log = True):
    if log:
        x = -torch.log(positives(Attrs).flatten())
    else:
        x = 1 - positives(Attrs).flatten().to(torch.float16)
    return x.sum()
def negative_loss(Attrs):
    """For two attributes
    """
    neg_pair = cos_sim(Attrs[0], Attrs[1])
    return torch.abs(neg_pair.flatten()).sum()
def train(attns, Z, optimizer, gamma = 2):
    attns.train()
    attrs = attns(Z)
    pos_loss = positive_loss(attrs)
    neg_loss = negative_loss(attrs)
    loss =  pos_loss + neg_loss * gamma
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()
    return pos_loss, neg_loss, loss

if __name__ == '__main__':
    saved = torch.load('all.pt')
    seed = 123
    ref = pd.read_csv('ref.csv', index_col = 0)
    Z = [z.hessian for z in saved]
    task2vec_dim = Z[0].shape[0]
    latent_dim = 128
    M = len(Z) # No of Task vectors
    Z_tensor = torch.tensor(Z)

    N = 2 # No of attributes
    attns = AttentionRaw((N, task2vec_dim)) #AttentionRawList(task2vec_dim, N)
    n_epochs = 4000
    gamma = 1.
    lr = 1e-1
    optimizer = torch.optim.Adam(attns.parameters(), lr = lr)
    for epoch in range(n_epochs):
        pos_loss, neg_loss, total_loss = train(attns, Z_tensor, optimizer, gamma = gamma )
        if epoch%200 == 0:
            print(f"Epoch: {epoch}, pos_loss: {pos_loss}, neg_loss: {neg_loss}, total_loss: {total_loss}")

    attrs = attns(Z_tensor)
    task_attributes = attrs.transpose(0,1)
    task_attribute_magnitudes = task_attributes.sum(-1)
    data = pd.DataFrame(task_attribute_magnitudes,  columns =['x','y'])
    data['lang'] = ref['lang']
    data['attr'] = ref['attr']
    data['name'] = ref['name']
    import seaborn as sns
    sns.set(rc={'figure.figsize':(15,15)})
    sns_plot = sns.scatterplot(data=data, x="x", y="y", hue = "attr", style="lang", s = 300)
    t = datetime.datetime.now()
    name = t.strftime('%I_%M%p_%B%d')
    sns_plot.get_figure().savefig(f"viz/soft-{n_epochs}-{lr}-{name}.png", dpi = 400)


text:
#!/usr/bin/python3
# -*- coding: utf-8 -*-
"""
OLIVIA MFSK Modulator
@author: Federico Santandrea

Source code is lightly commented for clarity.
Please read the document to better understand the various processes.
"""

# Libraries
import os, sys
import numpy as np
import sounddevice as sd
from numpy.fft import fft

# Hardcoded parameters
DEFAULT_PARAMS = "32/1000@1500"
SAMPLE_RATE = 8000
BLOCK_THRESHOLD = 24

# Global parameters
'''
Input:
- CENTER_FREQUENCY, SYMBOLS (# of tones), BANDWIDTH
Computed:
- spb: number of bits for a symbol
- fsep: frequency separation between tones, in Hz (also = baud rate)
- wlen: time separation between tones, in samples
'''

# Global objects
'''
- buf: audio sample input buffer
- sin: sounddevice InputStream for sample acquisition
'''

def __main__():
    '''
    Main program flow.
    '''
    welcomeAndSetup()
    initSound()

    syms = []

    while True:
        # Fetch new samples
        updateBuffer()

        sym = detectSymbol()
        syms.append(sym)

        if len(syms) == 64:
            # Enough symbols to decode a block
            if decodeAndPrintBlock(syms):
                # Block decoded successfully, waiting for a new one
                syms = []
            else:
                # Probably not a complete block, try rolling
                syms = syms[1:]

#

def welcomeAndSetup():
    '''
    Decodes command line parameters and prints a welcome message.
    '''
    global CENTER_FREQUENCY, SYMBOLS, BANDWIDTH
    global spb, fsep, wlen, slen

    if len(sys.argv) == 1:
        params = DEFAULT_PARAMS
    elif len(sys.argv) == 2:
        params = sys.argv[1]
    else:
        printUsageAndQuit()

    try:
        CENTER_FREQUENCY = int(params.split("@")[1])
        SYMBOLS = int(params.split("/")[0])
        BANDWIDTH = int(params.split("@")[0].split("/")[1])
        spb = int(np.log2(SYMBOLS))
        fsep = BANDWIDTH/SYMBOLS # = baud
        wlen = int(np.ceil(SAMPLE_RATE/fsep))

    except:
        printUsageAndQuit()

    if os.isatty(0):
        print("*** Olivia demodulator ***")
        print("(C) 2020 Federico Santandrea")
        print("federico.santandrea@studio.unibo.it")
        print()
        print("Starting Olivia demodulator at center " + str(CENTER_FREQUENCY) +
              "Hz, using " + str(SYMBOLS) + " tones over " + str(BANDWIDTH) + "Hz")
        print()

def printUsageAndQuit():
    '''
    Prints usage help if needed.
    '''
    print("usage: " + sys.argv[0] + " [syms/bandwidth@centerfrequency]")
    print("Example (default): " + sys.argv[0] + " " + DEFAULT_PARAMS)
    quit()

def initSound():
    '''
    Prepares global InputStream for sample acquisition.
    '''
    global sin, buf

    sin = sd.InputStream(samplerate=SAMPLE_RATE, blocksize=wlen,
        #channels=2, device=2, # stereo mix
        dtype=np.float32)
    sin.start()
    buf= np.zeros(wlen)

def updateBuffer():
    '''
    Acquires a new wlen-ful of samples from audio device.
    '''
    global buf

    (samples, of) = sin.read(wlen)
    buf = samples[:,0] # consider only one channel

def detectSymbol():
    '''
    Applies Fourier transform to audio buffer to detect
    symbol corresponding to sampled tone.

    Returns
    -------
    int
        Most likely symbol number.
    '''
    spectrum = np.abs(fft(buf))
    ix = CENTER_FREQUENCY - BANDWIDTH/2 + fsep/2
    measures = np.zeros(SYMBOLS)

    for i in range(0, SYMBOLS):
        ix += fsep
        measures[i] = spectrum[int(ix * wlen / SAMPLE_RATE)]
    mix = np.argmax(measures)

    return degray(mix)

def decodeAndPrintBlock(syms):
    '''
    Decodes a full block of 64 symbols, then prints it
    to standard output.
    '''
    w = np.zeros((spb, 64))

    # key = 0xE257E6D0291574EC
    key = np.flip(np.array(
          [1, 1, 1, 0, 0, 0, 1, 0,
           0, 1, 0, 1, 0, 1, 1, 1,
           1, 1, 1, 0, 0, 1, 1, 0,
           1, 1, 0, 1, 0, 0, 0, 0,
           0, 0, 1, 0, 1, 0, 0, 1,
           0, 0, 0, 1, 0, 1, 0, 1,
           0, 1, 1, 1, 0, 1, 0, 0,
           1, 1, 1, 0, 1, 1, 0, 0]))

    output = ""
    doubt = 0
    for i in range(0, spb):
        for j in range(0, 64):
            bit = (syms[j] >> ((i+j) % spb)) & 1
            if bit == 1:
                w[i,j] = -1
            else:
                w[i,j] = 1

        w[i,:] = w[i,:] * (-2*np.roll(key, -13*i)+1)
        w[i,:] = fwht(w[i,:])

        c = np.argmax(np.abs(w[i,:]))

        if abs(w[i,c]) < BLOCK_THRESHOLD:
            doubt += 1

        if w[i,c] < 0:
            c = c + 64
        if c != 0:
            output += chr(c)

    if doubt == 0:
        print(output, end="", flush=True)
        return True
    else:
        return False

#

def fwht(data):
    '''
    Fast Walsh-Hadamard transform.
    '''
    step = 1
    while step < len(data):
        for ptr in range(0, len(data), 2*step):
            for ptr2 in range(ptr, step+ptr):
                bit1 = data[ptr2]
                bit2 = data[ptr2+step]

                newbit1 = bit2
                newbit1 = newbit1 + bit1

                newbit2 = bit2
                newbit2 = newbit2 - bit1

                data[ptr2] = newbit1
                data[ptr2+step] = newbit2

        step *= 2
    return data

def degray(n):
    mask = n
    while mask != 0:
        mask >>= 1
        n ^= mask
    return n

#
#

__main__()

max_stars_repo_name: blikkentrommel/olivia-python

max_issues_repo_name: blikkentrommel/olivia-python

max_forks_repo_name: blikkentrommel/olivia-python

reason: import numpy,from numpy


text:
import numpy as np
import trimesh

def as_mesh(scene_or_mesh):
    """
    Convert a possible scene to a mesh.

    If conversion occurs, the returned mesh has only vertex and face data.
    """
    if isinstance(scene_or_mesh, trimesh.Scene):
        if len(scene_or_mesh.geometry) == 0:
            mesh = None  # empty scene
        else:
            # we lose texture information here
            mesh = trimesh.util.concatenate(
                tuple(trimesh.Trimesh(vertices=g.vertices, faces=g.faces)
                    for g in scene_or_mesh.geometry.values()))
    else:
        assert(isinstance(scene_or_mesh, trimesh.Trimesh))
        mesh = scene_or_mesh
    return mesh

def merge_meshes(mesh_list):
    vertices = []
    faces = []
    face_offset = 0
    for mesh in mesh_list:
        num_vertices = len(mesh.vertices)
        vertices.append(mesh.vertices)
        faces.append(mesh.faces + face_offset)
        face_offset += num_vertices
    vertices = np.concatenate(vertices, axis=0)
    faces = np.concatenate(faces, axis=0)
    merged_mesh = trimesh.Trimesh(vertices=vertices, faces=faces)
    return merged_mesh

def scale_to_unit_sphere(mesh):
    if isinstance(mesh, trimesh.Scene):
        mesh = mesh.dump().sum()

    vertices = mesh.vertices - mesh.bounding_box.centroid
    distances = np.linalg.norm(vertices, axis=1)
    vertices /= np.max(distances)

    return trimesh.Trimesh(vertices=vertices, faces=mesh.faces)


text:
# -*- coding: utf-8 -*-
from PIL import Image, ImageDraw, ImageFont
import numpy as np
import os
import random
import imageio
import glob
import sys

from deap import base
from deap import creator
from deap import tools
import deap_tools as extools

class Toolbox(base.Toolbox):
    def __init__(self, src_img_path, dst_alph, model):
        super(Toolbox, self).__init__()
        self.src_img_path = src_img_path
        self.dst_alph = dst_alph
        self.model = model
        self._set()

    def _load_img_as_np(self, path):
        img_pil = Image.open(path)
        if img_pil.size != (200, 200):
            img_pil = img_pil.resize((200, 200))
        img_np = np.asarray(img_pil)
        img_np.flags.writeable = True
        return img_np

    def eval_char(self, individual):
        pred = self.model.get_likelihoods(img_np=individual)
        dst_alph_num = ord(self.dst_alph) - 65
        return pred[dst_alph_num],

    def _set(self):
        creator.create('FitnessMax', base.Fitness, weights=(1.0,))
        creator.create('Individual', np.ndarray, fitness=creator.FitnessMax)
        self.register('attr_img', self._load_img_as_np, self.src_img_path)
        self.register('individual', tools.initIterate, creator.Individual, self.attr_img)
        self.register('population', tools.initRepeat, list, self.individual)
        self.register('evaluate', self.eval_char)
        self.register('mate', extools.cxTwoPointImg)
        self.register('mutate', extools.mutFlipBitImg, indpb=0.05)
        self.register('select', tools.selTournament, tournsize=3)

class AdversarialCharacter():
    def __init__(self, src_img_path, src_alph, dst_alph, dst_path,
                 cxpb, mutpb, ngen, npop, breakacc, model):
        self.src_alph = src_alph
        self.dst_alph = dst_alph
        self.dst_root_path = dst_path
        self.dst_best_path = os.path.join(self.dst_root_path, 'best')
        self.cxpb = cxpb
        self.mutpb = mutpb
        self.ngen = ngen
        self.npop = npop
        self.breakacc = breakacc
        self.accuracies = []
        self.model = model
        self.toolbox = Toolbox(src_img_path=src_img_path, dst_alph=dst_alph, model=model)
        self._make_dst_dir()

    def _make_dst_dir(self):
        if not os.path.exists(self.dst_root_path):
            os.mkdir(self.dst_root_path)
        if not os.path.exists(self.dst_best_path):
            os.mkdir(self.dst_best_path)

    def _save_img(self, filename, img_np):
        img_pil = Image.fromarray(np.uint8(img_np))
        img_pil.save(os.path.join(self.dst_best_path, filename), 'PNG')

    def _log_accuracies(self, best_ind_np):
        src_alph_score, dst_alph_score = self.model.eval_char_with_before(best_ind_np, self.src_alph, self.dst_alph)
        self.accuracies.append((src_alph_score, dst_alph_score))

    def _os_font_path(self):
        if sys.platform == 'darwin':
            return '/Library/Fonts/Arial.ttf'
        if sys.platform == 'linux':
            return '/usr/share/fonts/truetype/msttcorefonts/Arial.ttf'
        if sys.platform == 'win32':
            return 'C:\Windows\Fonts\arial.ttf'

    def train(self):
        # 初期集団を生成
        self.pop = self.toolbox.population(n=self.npop)

        # 初期集団の個体を評価
        fitnesses = list(map(self.toolbox.evaluate, self.pop))
        for ind, fit in zip(self.pop, fitnesses):
            ind.fitness.values = fit

        for g in range(self.ngen):
            print ('{0}世代'.format(g + 1))

            # 選択
            offspring = self.toolbox.select(self.pop, len(self.pop))
            offspring = list(map(self.toolbox.clone, offspring))

            # 交叉
            for child1, child2 in zip(offspring[::2], offspring[1::2]):
                if random.random() < self.cxpb:
                    self.toolbox.mate(child1, child2)
                    del child1.fitness.values
                    del child2.fitness.values

            # 変異
            for mutant in offspring:
                if random.random() < self.mutpb:
                    self.toolbox.mutate(mutant)
                    del mutant.fitness.values

            # 適合度が計算されていない個体を集めて適合度を計算
            invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
            fitnesses = map(self.toolbox.evaluate, invalid_ind)
            for ind, fit in zip(invalid_ind, fitnesses):
                ind.fitness.values = fit

            # 次世代群をoffspringにする
            self.pop[:] = offspring

            # すべての適合度を配列にする
            fits = [ind.fitness.values[0] for ind in self.pop]

            print ('Max: {0:013.10f} %'.format(max(fits) * 100))
            print ('Min: {0:013.10f} %'.format(min(fits) * 100))

            best_ind_np = tools.selBest(self.pop, 1)[0]
            self._save_img(str(g) + '.png', best_ind_np)
            self._log_accuracies(best_ind_np)

            self.finish_g = g + 1
            if max(fits) >= self.breakacc:
                break

    def save_animation(self, is_acc=True):
        imgs = []
        for i in range(self.finish_g):
            path = os.path.join(self.dst_best_path, '{0}.png'.format(i))
            img_char = Image.open(path)
            if is_acc:
                img = Image.new('RGB', (img_char.size[0], img_char.size[1] + 14), (255, 255, 255))
                draw = ImageDraw.Draw(img)
                font = ImageFont.truetype(self._os_font_path(), 14)
                src_acc_text = '{0}:{1:06.2f}%'.format(self.src_alph, self.accuracies[i][0] * 100)
                dst_acc_text = '{0}:{1:06.2f}%'.format(self.dst_alph, self.accuracies[i][1] * 100)
                iter_text = '{0:03d}/{1:03d}'.format(i + 1, self.finish_g)
                draw.text((0, img_char.size[1]), iter_text, font=font, fill='#000000')
                draw.text((58, img_char.size[1]), src_acc_text, font=font, fill='#000000')
                draw.text((130, img_char.size[1]), dst_acc_text, font=font, fill='#000000')
                img.paste(img_char, (0, 0))
            img_np = np.array(img)
            imgs.append(img_np)
        imageio.mimsave(os.path.join(self.dst_root_path, 'output.gif'), imgs, duration=0.5)
        print('save gif animation')

    def save_log(self):
        with open(os.path.join(self.dst_root_path, 'log.csv'), 'w') as log_file:
            log_file.write('{0},{1}\n'.format(self.src_alph, self.dst_alph))
            for i in range(self.finish_g):
                log_file.write('{0:012.10f},{1:012.10f}\n'.format(self.accuracies[i][0], self.accuracies[i][1]))
        print('save log')


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

# Max-Planck-Gesellschaft zur Förderung der Wissenschaften e.V. (MPG) is
# holder of all proprietary rights on this computer program.
# You can only use this computer program if you have closed
# a license agreement with MPG or you get the right to use the computer
# program from someone who is authorized to grant you that right.
# Any use of the computer program without a valid license is prohibited and
# liable to prosecution.
#
# Copyright©2019 Max-Planck-Gesellschaft zur Förderung
# der Wissenschaften e.V. (MPG). acting on behalf of its Max Planck Institute
# for Intelligent Systems. All rights reserved.
#
# Contact: ps-license@tuebingen.mpg.de

import time
import torch
import shutil
import logging
import numpy as np
import os.path as osp
from progress.bar import Bar

from lib.core.config import BASE_DATA_DIR
from lib.utils.utils import move_dict_to_device, AverageMeter, check_data_pararell

from lib.utils.eval_utils import (
    compute_accel,
    compute_error_accel,
    compute_error_verts,
    batch_compute_similarity_transform_torch,
)

logger = logging.getLogger(__name__)

class Trainer():
    def __init__(
            self,
            data_loaders,
            generator,
            motion_discriminator,
            gen_optimizer,
            dis_motion_optimizer,
            dis_motion_update_steps,
            end_epoch,
            criterion,
            start_epoch=0,
            lr_scheduler=None,
            motion_lr_scheduler=None,
            device=None,
            writer=None,
            debug=False,
            debug_freq=1000,
            logdir='output',
            resume=None,
            performance_type='min',
            num_iters_per_epoch=1000,
    ):

        # exclude motion discriminator
        # self.train_2d_loader, self.train_3d_loader, self.disc_motion_loader, self.valid_loader = data_loaders
        self.train_2d_loader, self.train_3d_loader, self.valid_loader = data_loaders

        # exclude motion discriminator
        # self.disc_motion_iter = iter(self.disc_motion_loader)

        self.train_2d_iter = self.train_3d_iter = None

        if self.train_2d_loader:
            self.train_2d_iter = iter(self.train_2d_loader)

        if self.train_3d_loader:
            self.train_3d_iter = iter(self.train_3d_loader)

        # Models and optimizers
        self.generator = generator
        self.gen_optimizer = gen_optimizer

        self.motion_discriminator = motion_discriminator
        self.dis_motion_optimizer = dis_motion_optimizer

        # Training parameters
        self.start_epoch = start_epoch
        self.end_epoch = end_epoch
        self.criterion = criterion
        self.lr_scheduler = lr_scheduler
        self.motion_lr_scheduler = motion_lr_scheduler
        self.device = device
        self.writer = writer
        self.debug = debug
        self.debug_freq = debug_freq
        self.logdir = logdir

        self.dis_motion_update_steps = dis_motion_update_steps

        self.performance_type = performance_type
        self.train_global_step = 0
        self.valid_global_step = 0
        self.epoch = 0
        self.best_performance = float('inf') if performance_type == 'min' else -float('inf')

        self.evaluation_accumulators = dict.fromkeys(['pred_j3d', 'target_j3d', 'target_theta', 'pred_verts'])

        self.num_iters_per_epoch = num_iters_per_epoch

        if self.writer is None:
            from torch.utils.tensorboard import SummaryWriter
            self.writer = SummaryWriter(log_dir=self.logdir)

        if self.device is None:
            self.device = 'cuda' if torch.cuda.is_available() else 'cpu'

        # Resume from a pretrained model
        if resume is not None:
            self.resume_pretrained(resume)

    def train(self):
        # Single epoch training routine

        losses = AverageMeter()
        kp_2d_loss = AverageMeter()
        kp_3d_loss = AverageMeter()
        # accel_loss = AverageMeter()

        timer = {
            'data': 0,
            'forward': 0,
            'loss': 0,
            'backward': 0,
            'batch': 0,
        }

        self.generator.train()
        self.motion_discriminator.train()

        start = time.time()

        summary_string = ''

        bar = Bar(f'Epoch {self.epoch + 1}/{self.end_epoch}', fill='#', max=self.num_iters_per_epoch)

        for i in range(self.num_iters_per_epoch):
            # Dirty solution to reset an iterator
            target_2d = target_3d = None
            if self.train_2d_iter:
                try:
                    target_2d = next(self.train_2d_iter)
                except StopIteration:
                    self.train_2d_iter = iter(self.train_2d_loader)
                    target_2d = next(self.train_2d_iter)

                move_dict_to_device(target_2d, self.device)

            if self.train_3d_iter:
                try:
                    target_3d = next(self.train_3d_iter)
                except StopIteration:
                    self.train_3d_iter = iter(self.train_3d_loader)
                    target_3d = next(self.train_3d_iter)

                move_dict_to_device(target_3d, self.device)

            # exclude motion discriminator
            # real_body_samples = real_motion_samples = None
            # try:
            #     real_motion_samples = next(self.disc_motion_iter)
            # except StopIteration:
            #     self.disc_motion_iter = iter(self.disc_motion_loader)
            #     real_motion_samples = next(self.disc_motion_iter)
            # move_dict_to_device(real_motion_samples, self.device)

            # <======= Feedforward generator and discriminator
            if target_2d and target_3d:
                inp = torch.cat((target_2d['features'], target_3d['features']), dim=0).cuda()
            elif target_3d:
                inp = target_3d['features'].cuda()
            else:
                inp = target_2d['features'].cuda()

            timer['data'] = time.time() - start
            start = time.time()

            preds, scores = self.generator(inp, is_train=True)

            timer['forward'] = time.time() - start
            start = time.time()

            gen_loss, loss_dict = self.criterion(
                generator_outputs=preds,
                data_2d=target_2d,
                data_3d=target_3d,
                scores=scores
            )

            # exclude motion discriminator
            # gen_loss, motion_dis_loss, loss_dict = self.criterion(
            #     generator_outputs=preds,
            #     data_2d=target_2d,
            #     data_3d=target_3d,
            #     scores=scores,
            #     data_body_mosh=real_body_samples,
            #     data_motion_mosh=real_motion_samples,
            #     motion_discriminator=self.motion_discriminator,
            # )
            # =======>

            timer['loss'] = time.time() - start
            start = time.time()

            # <======= Backprop generator and discriminator
            self.gen_optimizer.zero_grad()
            gen_loss.backward()
            self.gen_optimizer.step()

            # exclude motion discriminator
            # if self.train_global_step % self.dis_motion_update_steps == 0:
            #     self.dis_motion_optimizer.zero_grad()
            #     motion_dis_loss.backward()
            #     self.dis_motion_optimizer.step()
            # =======>

            # <======= Log training info
            total_loss = gen_loss
            # exclude motion discriminator
            # total_loss = gen_loss + motion_dis_loss

            losses.update(total_loss.item(), inp.size(0))
            kp_2d_loss.update(loss_dict['loss_kp_2d'].item(), inp.size(0))
            kp_3d_loss.update(loss_dict['loss_kp_3d'].item(), inp.size(0))

            timer['backward'] = time.time() - start
            timer['batch'] = timer['data'] + timer['forward'] + timer['loss'] + timer['backward']
            start = time.time()

            summary_string = f'({i + 1}/{self.num_iters_per_epoch}) | Total: {bar.elapsed_td} | ' \
                             f'ETA: {bar.eta_td:} | loss: {losses.avg:.2f} | 2d: {kp_2d_loss.avg:.2f} ' \
                             f'| 3d: {kp_3d_loss.avg:.2f} '

            for k, v in loss_dict.items():
                summary_string += f' | {k}: {v:.3f}'
                self.writer.add_scalar('train_loss/'+k, v, global_step=self.train_global_step)

            for k,v in timer.items():
                summary_string += f' | {k}: {v:.2f}'

            self.writer.add_scalar('train_loss/loss', total_loss.item(), global_step=self.train_global_step)

            if self.debug:
                print('==== Visualize ====')
                from lib.utils.vis import batch_visualize_vid_preds
                video = target_3d['video']
                dataset = 'spin'
                vid_tensor = batch_visualize_vid_preds(video, preds[-1], target_3d.copy(),
                                                       vis_hmr=False, dataset=dataset)
                self.writer.add_video('train-video', vid_tensor, global_step=self.train_global_step, fps=10)

            self.train_global_step += 1
            bar.suffix = summary_string
            bar.next()

            if torch.isnan(total_loss):
                exit('Nan value in loss, exiting!...')
            # =======>

        bar.finish()

        logger.info(summary_string)

    def validate(self):
        self.generator.eval()

        start = time.time()

        summary_string = ''

        bar = Bar('Validation', fill='#', max=len(self.valid_loader))

        if self.evaluation_accumulators is not None:
            for k,v in self.evaluation_accumulators.items():
                self.evaluation_accumulators[k] = []

        J_regressor = torch.from_numpy(np.load(osp.join(BASE_DATA_DIR, 'J_regressor_h36m.npy'))).float()

        for i, target in enumerate(self.valid_loader):

            move_dict_to_device(target, self.device)

            # <=============
            with torch.no_grad():
                inp = target['features']
                batch = len(inp)

                preds, _ = self.generator(inp, J_regressor=J_regressor)

                # convert to 14 keypoint format for evaluation
                n_kp = preds[-1]['kp_3d'].shape[-2]
                pred_j3d = preds[-1]['kp_3d'].view(-1, n_kp, 3).cpu().numpy()
                target_j3d = target['kp_3d'].view(-1, n_kp, 3).cpu().numpy()
                pred_verts = preds[-1]['verts'].view(-1, 6890, 3).cpu().numpy()
                target_theta = target['theta'].view(-1, 85).cpu().numpy()

                self.evaluation_accumulators['pred_verts'].append(pred_verts)
                self.evaluation_accumulators['target_theta'].append(target_theta)

                self.evaluation_accumulators['pred_j3d'].append(pred_j3d)
                self.evaluation_accumulators['target_j3d'].append(target_j3d)
            # =============>

            # <============= DEBUG
            if self.debug and self.valid_global_step % self.debug_freq == 0:
                from lib.utils.vis import batch_visualize_vid_preds
                video = target['video']
                dataset = 'common'
                vid_tensor = batch_visualize_vid_preds(video, preds[-1], target, vis_hmr=False, dataset=dataset)
                self.writer.add_video('valid-video', vid_tensor, global_step=self.valid_global_step, fps=10)
            # =============>

            batch_time = time.time() - start

            summary_string = f'({i + 1}/{len(self.valid_loader)}) | batch: {batch_time * 10.0:.4}ms | ' \
                             f'Total: {bar.elapsed_td} | ETA: {bar.eta_td:}'

            self.valid_global_step += 1
            bar.suffix = summary_string
            bar.next()

        bar.finish()

        logger.info(summary_string)

    def fit(self):

        for epoch in range(self.start_epoch, self.end_epoch):
            self.epoch = epoch
            self.train()
            self.validate()
            performance = self.evaluate()

            if self.lr_scheduler is not None:
                self.lr_scheduler.step(performance)

            if self.motion_lr_scheduler is not None:
                self.motion_lr_scheduler.step(performance)

            # log the learning rate
            for param_group in self.gen_optimizer.param_groups:
                print(f'Learning rate {param_group["lr"]}')
                self.writer.add_scalar('lr/gen_lr', param_group['lr'], global_step=self.epoch)

            for param_group in self.dis_motion_optimizer.param_groups:
                print(f'Learning rate {param_group["lr"]}')
                self.writer.add_scalar('lr/dis_lr', param_group['lr'], global_step=self.epoch)

            logger.info(f'Epoch {epoch+1} performance: {performance:.4f}')

            self.save_model(performance, epoch)

        self.writer.close()

    def save_model(self, performance, epoch):
        save_dict = {
            'epoch': epoch,
            'gen_state_dict': self.generator.state_dict(),
            'performance': performance,
            'gen_optimizer': self.gen_optimizer.state_dict(),
            'disc_motion_state_dict': self.motion_discriminator.state_dict(),
            'disc_motion_optimizer': self.dis_motion_optimizer.state_dict(),
        }

        filename = osp.join(self.logdir, 'checkpoint.pth.tar')
        torch.save(save_dict, filename)

        if self.performance_type == 'min':
            is_best = performance < self.best_performance
        else:
            is_best = performance > self.best_performance

        if is_best:
            logger.info('Best performance achived, saving it!')
            self.best_performance = performance
            shutil.copyfile(filename, osp.join(self.logdir, 'model_best.pth.tar'))

            with open(osp.join(self.logdir, 'best.txt'), 'w') as f:
                f.write(str(float(performance)))

    def resume_pretrained(self, model_path):
        if osp.isfile(model_path):
            checkpoint = torch.load(model_path)
            self.start_epoch = checkpoint['epoch']
            self.generator.load_state_dict(checkpoint['gen_state_dict'])
            self.gen_optimizer.load_state_dict(checkpoint['gen_optimizer'])
            self.best_performance = checkpoint['performance']

            if 'disc_motion_optimizer' in checkpoint.keys():
                self.motion_discriminator.load_state_dict(checkpoint['disc_motion_state_dict'])
                self.dis_motion_optimizer.load_state_dict(checkpoint['disc_motion_optimizer'])

            logger.info(f"=> loaded checkpoint '{model_path}' "
                  f"(epoch {self.start_epoch}, performance {self.best_performance})")
        else:
            logger.info(f"=> no checkpoint found at '{model_path}'")

    def evaluate(self):

        for k, v in self.evaluation_accumulators.items():
            self.evaluation_accumulators[k] = np.vstack(v)

        pred_j3ds = self.evaluation_accumulators['pred_j3d']
        target_j3ds = self.evaluation_accumulators['target_j3d']

        pred_j3ds = torch.from_numpy(pred_j3ds).float()
        target_j3ds = torch.from_numpy(target_j3ds).float()

        print(f'Evaluating on {pred_j3ds.shape[0]} number of poses...')
        pred_pelvis = (pred_j3ds[:,[2],:] + pred_j3ds[:,[3],:]) / 2.0
        target_pelvis = (target_j3ds[:,[2],:] + target_j3ds[:,[3],:]) / 2.0

        pred_j3ds -= pred_pelvis
        target_j3ds -= target_pelvis

        errors = torch.sqrt(((pred_j3ds - target_j3ds) ** 2).sum(dim=-1)).mean(dim=-1).cpu().numpy()
        S1_hat = batch_compute_similarity_transform_torch(pred_j3ds, target_j3ds)
        errors_pa = torch.sqrt(((S1_hat - target_j3ds) ** 2).sum(dim=-1)).mean(dim=-1).cpu().numpy()
        pred_verts = self.evaluation_accumulators['pred_verts']
        target_theta = self.evaluation_accumulators['target_theta']

        m2mm = 1000

        pve = np.mean(compute_error_verts(target_theta=target_theta, pred_verts=pred_verts)) * m2mm
        accel = np.mean(compute_accel(pred_j3ds)) * m2mm
        accel_err = np.mean(compute_error_accel(joints_pred=pred_j3ds, joints_gt=target_j3ds)) * m2mm
        mpjpe = np.mean(errors) * m2mm
        pa_mpjpe = np.mean(errors_pa) * m2mm

        eval_dict = {
            'mpjpe': mpjpe,
            'pa-mpjpe': pa_mpjpe,
            'accel': accel,
            'pve': pve,
            'accel_err': accel_err
        }

        log_str = f'Epoch {self.epoch}, '
        log_str += ' '.join([f'{k.upper()}: {v:.4f},'for k,v in eval_dict.items()])
        logger.info(log_str)

        for k,v in eval_dict.items():
            self.writer.add_scalar(f'error/{k}', v, global_step=self.epoch)

        # return accel_err
        return pa_mpjpe


text:
"""Runs the environments located in flow/benchmarks.
The environment file can be modified in the imports to change the environment
this runner script is executed on. This file runs the DDPG algorithm in rllib
and utilizes the hyper-parameters specified in:
Proximal Policy Optimization Algorithms by Schulman et. al.
"""
import json
import argparse
from itertools import product
import numpy as np

import ray
try:
    from ray.rllib.agents.agent import get_agent_class
except ImportError:
    from ray.rllib.agents.registry import get_agent_class
from ray.tune import Experiment, run_experiments
from ray.tune.registry import register_env

from flow.utils.registry import make_create_env
from flow.utils.rllib import FlowParamsEncoder

EXAMPLE_USAGE = """
example usage:
    python ddpg_runner.py grid0
Here the arguments are:
benchmark_name - name of the benchmark to run
num_cpus - number of cpus to use for training
"""

parser = argparse.ArgumentParser(
    formatter_class=argparse.RawDescriptionHelpFormatter,
    description="[Flow] Evaluates a Flow Garden solution on a benchmark.",
    epilog=EXAMPLE_USAGE)

# required input parameters
parser.add_argument(
    "--benchmark_name", type=str, help="File path to solution environment.")

parser.add_argument(
    "--exp_tag", type=str, help="experiment tag")

# optional input parameters
parser.add_argument(
    '--num_cpus',
    type=int,
    default=63,
    help="The number of cpus to use.")

def on_episode_start(info):
    episode = info["episode"]
    episode.user_data["cost1"] = []
    episode.user_data["cost2"] = []
    episode.user_data["mean_vel"] = []
    episode.user_data["outflow"] = []

def on_episode_step(info):
    episode = info["episode"]
    agent_ids = episode.agent_rewards.keys()
    infos = [episode.last_info_for(id_[0]) for id_ in agent_ids]
    cost1, cost2, mean_vel, outflow = 0, 0, 0, 0
    if len(infos) != 0:
        cost1 = np.mean([info['cost1'] for info in infos])
        cost2 = np.mean([info['cost2'] for info in infos])
        mean_vel = np.mean([info['mean_vel'] for info in infos])
        outflow = np.mean([info['outflow'] for info in infos])
    episode.user_data["cost1"].append(cost1)
    episode.user_data["cost2"].append(cost2)
    episode.user_data["mean_vel"].append(mean_vel)
    episode.user_data["outflow"].append(outflow)

def on_episode_end(info):
    episode = info["episode"]
    cost1 = np.sum(episode.user_data["cost1"])
    cost2 = np.sum(episode.user_data["cost2"])
    mean_vel = np.mean(episode.user_data["mean_vel"])
    outflow = np.mean(episode.user_data["outflow"][-500:])  # 1/3 of the whole steps
    episode.custom_metrics["cost1"] = cost1
    episode.custom_metrics["cost2"] = cost2
    episode.custom_metrics["system_level_velocity"] = mean_vel
    episode.custom_metrics["outflow_rate"] = outflow

if __name__ == "__main__":
    args = parser.parse_args()
    benchmark_name = args.benchmark_name
    # number of parallel workers
    num_cpus = args.num_cpus

    # Import the benchmark and fetch its flow_params
    benchmark = __import__(
        "flow.benchmarks.%s" % benchmark_name, fromlist=["flow_params"])
    flow_params = benchmark.flow_params

    # initialize a ray instance
    ray.init()

    alg_run = "DDPG"

    horizon = flow_params["env"].horizon
    sim_step = flow_params["sim"].sim_step

    agent_cls = get_agent_class(alg_run)
    # use almost defalt config
    config = agent_cls._default_config.copy()
    config["num_workers"] = num_cpus
    config["train_batch_size"] = 3000
    config["sample_batch_size"] = 750
    config["horizon"] = horizon
    config['timesteps_per_iteration'] = int(horizon / sim_step)
    config['clip_actions'] = False  # FIXME(ev) temporary ray bug
    config["observation_filter"] = "NoFilter"
    config["batch_mode"] = "complete_episodes"

    # save the flow params for replay
    flow_json = json.dumps(
        flow_params, cls=FlowParamsEncoder, sort_keys=True, indent=4)
    config['env_config']['flow_params'] = flow_json
    config['env_config']['run'] = alg_run

    config['callbacks']['on_episode_start'] = ray.tune.function(on_episode_start)
    config['callbacks']['on_episode_step'] = ray.tune.function(on_episode_step)
    config['callbacks']['on_episode_end'] = ray.tune.function(on_episode_end)

    # tunning parameters
    e2_list = [0.1, 0.3]
    e3_list = [0.0]
    t_min = [10.0]

    env_name_list = []
    i = 0
    for e2, e3, t in product(e2_list, e3_list, t_min):
        flow_params["env"].additional_params["eta1"] = 1.0# e[0]
        flow_params["env"].additional_params["eta2"] = e2 # e[1]
        flow_params["env"].additional_params["eta3"] = e3 # e[2]
        # flow_params["env"].additional_params["reward_scale"] = rew
        flow_params["env"].additional_params["t_min"] = t
        flow_params["env"].additional_params["buf_length"] = 1

        # get the env name and a creator for the environment
        create_env, env_name = make_create_env(params=flow_params, version=0)
        env_name = env_name + '_[eta1, eta2, eta3]:[{}, {}, {}]'.format(1.0, e2, e3) + '_t_min:{}'.format(t)
        env_name_list.append(env_name)
        # Register as rllib env
        register_env(env_name, create_env)

    exp_list = []
    for env_name in env_name_list:
        exp_tag = {
            "run": alg_run,
            "env": env_name,
            "config": {
                **config
            },
            "checkpoint_freq": 25,
            "max_failures": 999,
            "stop": {
                "training_iteration": 50
            },
            "num_samples": 1,
        }
        exp_list.append(Experiment.from_json(args.exp_tag, exp_tag))

    trials = run_experiments(
        experiments=exp_list
    )