text:
#!/usr/bin/env python

import os
import json
import pickle
import numpy as np
import pandas as pd
from spacy.lang import en

from sklearn.manifold import TSNE
from sklearn.decomposition import PCA

from bokeh.palettes import brewer
from bokeh.plotting import figure, show, output_file, ColumnDataSource
from bokeh.models import HoverTool, WheelZoomTool, PanTool, BoxZoomTool, ResetTool, value

def get_pca_vectors(input_vec, n_components=200):
    pca = PCA(n_components=n_components)
    return pca.fit_transform(input_vec)

def get_tsne_vectors(input_vec, perplexity=40, early_exaggeration=10,
                n_iter=2500, v=2, serialize=False, export_dir=None, name=''):
    tsne = TSNE(perplexity=perplexity, early_exaggeration=early_exaggeration,
        n_iter=n_iter, verbose=v)
    tsne_vectors = tsne.fit_transform(input_vec)
    # save t-SNE model and vectors to disk
    if serialize:
        tsne_model_path = os.path.join(export_dir, 'tsne_' + name + '_model.pkl')
        with open(tsne_model_path, 'wb') as out:
            pickle.dump(tsne, out)
        tsne_vecs_path = os.path.join(export_dir, 'tsne_' + name + '_vectors.npy')
        np.save(tsne_vecs_path, tsne_vectors)
        print('t-SNE model and vectors saved at', os.path.realpath(export_dir))
    return tsne_vectors

def filter_unspecified(bvecs, tvecs, post_tags):
    new_bvecs = None
    new_tvecs = None
    new_post_tags = []
    for idx, tag in enumerate(post_tags):
        if tag != 'unspecified':
            new_post_tags.append(tag)
            if new_bvecs is None:
                new_bvecs = bvecs[idx]
                new_tvecs = tvecs[idx]
            else:
                new_bvecs = np.vstack((new_bvecs, bvecs[idx]))
                new_tvecs = np.vstack((new_tvecs, tvecs[idx]))
    return new_bvecs, new_tvecs, new_post_tags

def sort_by_score(data, metadata_path):
    with open(metadata_path, 'rb') as _in:
        metadata = json.load(_in)

    scores = []
    for post in metadata:
        scores.append(post['Score'])
    sorted_idx = np.argsort(scores)

    sorted_data = None
    if isinstance(data, np.ndarray):
        sorted_data = np.zeros(data.shape, dtype=data.dtype)
        for ii, idx in enumerate(sorted_idx[::-1]):
            sorted_data[ii] = data[idx]
    elif isinstance(data, list):
        sorted_data = []
        for idx in sorted_idx[::-1]:
            sorted_data.append(data[idx])

    return sorted_data

def assign_colors(post_tags, highlight_tags):
    palette = brewer['Set1'][len(highlight_tags)]
    colors = []
    tag_names = []
    for tag in post_tags:
        if tag in highlight_tags:
            tag_names.append(tag)
            colors.append(palette[highlight_tags.index(tag)])
        else:
            tag_names.append(tag)
            colors.append('#d3d3d3')
    return colors, tag_names

def build_legend(post_tags, highlight_tags):
    legend = []
    for tag in post_tags:
        if tag in highlight_tags:
            legend.append(tag)
        else:
            legend.append('unspecified')
    return legend

def posts_with_specific_tags(metadata_path, highlight_tags):
    with open(metadata_path, 'rb') as _in:
        metadata = json.load(_in)
    post_tags = []
    for post in metadata:
        no_tag = True
        for tag in post['Tags']:
            if tag in highlight_tags and len(post['Tags']) <= 5:
                no_tag = False
                post_tags.append(tag)
                break
        if no_tag:
            #post_tags.append('unspecified')
            post_tags.append(', '.join(post['Tags']))
    return post_tags

def bokeh_plot(dataframe, export_dir, name, plot_size=850):
    png_path = os.path.join(export_dir, 't-SNE_' + name + '_vectors.png')
    html_path = os.path.join(export_dir, 't-SNE_' + name + '_vectors.html')
    title = '[t-SNE] Java/StackOverflow ' + name + ' vectors in 2D space'
    output_file(html_path, title=title)
    # plot tools
    hover = HoverTool(tooltips = '@tag')
    tools = [hover, WheelZoomTool(), BoxZoomTool(), ResetTool(), PanTool()]
    #tools = [WheelZoomTool(), BoxZoomTool(), ResetTool(), PanTool()]
    # add DataFrame as a ColumnDataSource for Bokeh
    plot_data = ColumnDataSource(dataframe)
    # create the plot and configure the title, dimensions, and tools
    tsne_plot = figure(
                        title=title,
                        tools=tools,
                        plot_width = plot_size,
                        plot_height = plot_size)
    tsne_plot.circle(
                    'x_coord', 'y_coord',
                    source=plot_data, color='color', legend='legend',
                    line_alpha=0.2, fill_alpha=0.5, size=10, hover_line_color='black')
    # configure visual elements of the plot
    tsne_plot.title.text_font_size = value('16pt')
    tsne_plot.xaxis.visible = False
    tsne_plot.yaxis.visible = False
    tsne_plot.grid.grid_line_color = None
    tsne_plot.outline_line_color = None
    # engage!
    show(tsne_plot)

def main(vector_index_path, metadata_path, highlight_tags, export_dir,
                                        sample_size=4000, pca_smooth=True):
    with open(vector_index_path, 'rb') as _in:
        vec_index = pickle.load(_in)
    bvecs = vec_index['BodyVectors']
    tvecs = vec_index['TitleVectors']

    if pca_smooth:
        bvecs = get_pca_vectors(bvecs)
        tvecs = get_pca_vectors(tvecs)

    post_tags = posts_with_specific_tags(metadata_path, highlight_tags)

    # sort by score and get top sample_size rows
    bvecs = sort_by_score(bvecs, metadata_path)[:sample_size]
    tvecs = sort_by_score(tvecs, metadata_path)[:sample_size]
    post_tags = sort_by_score(post_tags, metadata_path)[:sample_size]
    #bvecs, tvecs, post_tags = filter_unspecified(bvecs, tvecs, post_tags)
    colors, tag_names = assign_colors(post_tags, highlight_tags)
    legend = build_legend(post_tags, highlight_tags)

    tsne_bvecs = get_tsne_vectors(bvecs, serialize=True, export_dir=export_dir,
        name='body')
    tsne_tvecs = get_tsne_vectors(tvecs, serialize=True, export_dir=export_dir,
        name='title')

    '''
    tsne_vecs_path = os.path.join(export_dir, 'tsne_body_vectors.npy')
    tsne_bvecs = np.load(tsne_vecs_path)
    tsne_vecs_path = os.path.join(export_dir, 'tsne_title_vectors.npy')
    tsne_tvecs = np.load(tsne_vecs_path)
    '''

    bdframe = pd.DataFrame(tsne_bvecs, columns=['x_coord', 'y_coord'])
    bdframe['color'] = colors
    bdframe['tag'] = tag_names
    bdframe['legend'] = legend
    bokeh_plot(bdframe, export_dir, 'body')

    tdframe = pd.DataFrame(tsne_tvecs, columns=['x_coord', 'y_coord'])
    tdframe['color'] = colors
    tdframe['tag'] = tag_names
    tdframe['legend'] = legend
    bokeh_plot(tdframe, export_dir, 'title')

if __name__ == '__main__':
    demo_tags = ['junit', 'jpa', 'jackson', 'multithreading', 'date']
    main('../../src/wordvec_model/index/ft_v0.1.2_post_index.pkl',
            '../../src/wordvec_model/index/metadata.json', demo_tags, '.')

max_stars_repo_name: nikosoik/stacksearch

max_issues_repo_name: nikosoik/stacksearch

max_forks_repo_name: nikosoik/stacksearch


text:
# Copyright (c) 2019 The Regents of the University of Michigan
# All rights reserved.
# This software is licensed under the BSD 3-Clause License.
"""Interpolate between pairs of quaternions.

The rowan package provides a simple interface to slerp, the standard method of
quaternion interpolation for two quaternions.
"""

import numpy as np

from ..functions import _validate_unit, conjugate, log, multiply, power

__all__ = ["slerp", "slerp_prime", "squad"]

def slerp(q0, q1, t, ensure_shortest=True):
    r"""Spherical linear interpolation between p and q.

    The `slerp formula <https://en.wikipedia.org/wiki/Slerp#Quaternion_Slerp>`_
    can be easily expressed in terms of the quaternion exponential (see
    :func:`rowan.exp`).

    Args:
        q0 ((..., 4) :class:`numpy.ndarray`):
            First array of quaternions.
        q1 ((..., 4) :class:`numpy.ndarray`):
            Second array of quaternions.
        t ((...) :class:`numpy.ndarray`):
            Interpolation parameter :math:`\in [0, 1]`
        ensure_shortest (bool):
            Flip quaternions to ensure we traverse the geodesic in the shorter
            (:math:`<180^{\circ}`) direction.

    .. note::

        Given inputs such that :math:`t\notin [0, 1]`, the values outside the
        range are simply assumed to be 0 or 1 (depending on which side of the
        interval they fall on).

    Returns:
        (..., 4) :class:`numpy.ndarray`: Interpolations between ``p`` and ``q``.

    Example::

        >>> import numpy as np
        >>> rowan.interpolate.slerp(
        ...     [[1, 0, 0, 0]], [[np.sqrt(2)/2, np.sqrt(2)/2, 0, 0]], 0.5)
        array([[0.92387953, 0.38268343, 0.        , 0.        ]])
    """
    _validate_unit(q0)
    _validate_unit(q1)
    t = np.clip(t, 0, 1)

    q0 = np.asarray(np.atleast_2d(q0))
    q1 = np.array(np.atleast_2d(q1))

    # Ensure that we turn the short way around
    if ensure_shortest:
        cos_theta = np.sum(q0 * q1, axis=-1)
        flip = cos_theta < 0
        q1[flip] *= -1

    return multiply(q0, power(multiply(conjugate(q0), q1), t))

def slerp_prime(q0, q1, t, ensure_shortest=True):
    r"""Compute the derivative of slerp.

    Args:
        q0 ((..., 4) :class:`numpy.ndarray`):
            First set of quaternions.
        q1 ((..., 4) :class:`numpy.ndarray`):
            Second set of quaternions.
        t ((...) :class:`numpy.ndarray`):
            Interpolation parameter :math:`\in [0, 1]`
        ensure_shortest (bool):
            Flip quaternions to ensure we traverse the geodesic in the shorter
            (:math:`<180^{\circ}`) direction

    Returns:
        (..., 4) :class:`numpy.ndarray`:
            The derivative of the interpolations between ``p`` and ``q``.

    Example::

        import numpy as np
        q_slerp_prime rowan.interpolate.slerp_prime(
            [[1, 0, 0, 0]], [[np.sqrt(2)/2, np.sqrt(2)/2, 0, 0]], 0.5)
    """
    _validate_unit(q0)
    _validate_unit(q1)
    t = np.clip(t, 0, 1)

    q0 = np.asarray(np.atleast_2d(q0))
    q1 = np.array(np.atleast_2d(q1))

    # Ensure that we turn the short way around
    if ensure_shortest:
        cos_theta = np.sum(q0 * q1, axis=-1)
        flip = cos_theta < 0
        q1[flip] *= -1

    return multiply(
        multiply(q0, power(multiply(conjugate(q0), q1), t)),
        log(multiply(conjugate(q0), q1)),
    )

def squad(p, a, b, q, t):
    r"""Cubically interpolate between p and q.

    The SQUAD formula is just a repeated application of Slerp between multiple
    quaternions as originally derived in [Shoemake85]_:

    .. math::
        \begin{equation}
            \textrm{squad}(p, a, b, q, t) = \textrm{slerp}(p, q, t)
            \left(\textrm{slerp}(p, q, t)^{-1}\textrm{slerp}(a, b, t)
            \right)^{2t(1-t)}
        \end{equation}

    .. [Shoemake85] Ken Shoemake. Animating rotation with quaternion curves.
        SIGGRAPH Comput. Graph., 19(3):245-254, July 1985.

    Args:
        p ((..., 4) :class:`numpy.ndarray`): First endpoint of interpolation.
        a ((..., 4) :class:`numpy.ndarray`): First control point of interpolation.
        b ((..., 4) :class:`numpy.ndarray`): Second control point of interpolation.
        q ((..., 4) :class:`numpy.ndarray`): Second endpoint of interpolation.
        t ((...) :class:`numpy.ndarray`): Interpolation parameter :math:`t \in [0, 1]`.

    Returns:
        (..., 4) :class:`numpy.ndarray`:
            Interpolations between ``p`` and ``q``.

    Example::

        >>> import numpy as np
        >>> rowan.interpolate.squad(
        ...     [1, 0, 0, 0], [np.sqrt(2)/2, np.sqrt(2)/2, 0, 0],
        ...     [0, np.sqrt(2)/2, np.sqrt(2)/2, 0],
        ...     [0, 0, np.sqrt(2)/2, np.sqrt(2)/2], 0.5)
        array([[0.64550177, 0.47254009, 0.52564058, 0.28937053]])
    """
    _validate_unit(p)
    _validate_unit(a)
    _validate_unit(b)
    _validate_unit(q)
    t = np.clip(t, 0, 1)

    return slerp(
        slerp(p, q, t, ensure_shortest=False),
        slerp(a, b, t, ensure_shortest=False),
        2 * t * (1 - t),
        ensure_shortest=False,
    )

ext: py

max_stars_repo_name: glotzerlab/rowan

max_issues_repo_name: glotzerlab/rowan

max_forks_repo_name: glotzerlab/rowan

reason: import numpy


text:
import numpy as np
import blarf
import os

natoms = 6

k = 6

ff = blarf.potential()

timemax = 10.0
stride = 100
timestep = 0.001

ff.init_ethylene(1.0,0.05,0.3)

r = np.zeros(natoms*3)

r[5] = 1.1
r[7] = 0.7
r[8] = -0.5
r[9] = 0.1
r[10] = -0.7
r[11] = -0.5
r[12] = 0.7
r[14] = 1.5
r[15] = -0.6
r[17] = 1.5

print "r ", r

ff.eval_pes(r)

print "e ", ff.get_energy()
print "f ", ff.get_forces()

vv = blarf.classical()

vv.set_numdims(natoms*3)
vv.set_positions(r)
vv.set_momenta(np.zeros(natoms*3))
vv.set_timestep(timestep)
vv.set_maxtime(timemax)
vv.set_num_rescale_steps(100)
vv.set_temperature(0.3)

data = blarf.dataset(natoms*3)

vv.propagate(ff,data,stride=stride)

print data.get_numpoints()
print data.get_positions()
print data.get_energies_exact()

thindata = data.thin_dataset(1)
thindata.set_internal_type("reciprical_bond")
thindata.compute_internals()

r2 = np.zeros(natoms*3)
r2[5] = 1.2
r2[6] = -0.4
r2[7] = 0.5
r2[8] = -0.6
r2[9] = 0.3
r2[10] = -0.5
r2[11] = -0.5
r2[12] = 0.7
r2[14] = 1.4
r2[15] = -0.6
r2[16] = 0.1
r2[17] = 1.6

vv2 = blarf.classical()

vv2.set_numdims(natoms*3)
vv2.set_positions(r2)
vv2.set_momenta(np.zeros(natoms*3))
vv2.set_timestep(timestep)
vv2.set_maxtime(timemax)
vv2.set_num_rescale_steps(100)
vv2.set_temperature(0.3)

data2 = blarf.dataset(natoms*3)

vv2.propagate(ff,data2,stride=stride)

thindata2 = data2.thin_dataset(1)
thindata2.set_internal_type("reciprical_bond")
thindata2.compute_internals()

#width_factors = [2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0]
#numtrials = 5
width_factors = [1.0]
numtrials = 1

if os.path.exists("mur.dat"):
    os.remove("mur.dat")

if os.path.exists("mur2.dat"):
    os.remove("mur2.dat")

iw = 0
for w in width_factors:
    murs = np.zeros(numtrials)
    murs2 = np.zeros(numtrials)
    for itrial in range(numtrials):

        clust = blarf.cluster(k,thindata.get_numinternals())

        clust.k_means_optimization(thindata)

        network = blarf.rbfn()

        network.set_width_factor(w)
        network.set_regularization_constant(1.0)
        #network.init_from_cluster_reciprical_bonds_traditionalrbf(clust)
        network.init_from_cluster_reciprical_bonds_onedimensional(clust)

        network.solve_weights(thindata)

        #filename = "network." + str(w) + "." + str(itrial) + ".hdf5"

        #network.h5_output(filename)

        murs[itrial] = thindata.get_mean_unsigned_residual()

        #thindata.h5_output("thindata.hdf5")

        network.compute_energies_approx(thindata2)

        murs2[itrial] = thindata2.get_mean_unsigned_residual()

        #thindata2.h5_output("thindata2.hdf5")

    f = open('mur2-eth.dat', 'a')

    s = str(w) + "  "
    for imur in range(murs.size):
        s = s + str(murs[imur]) + "  "
    s = s + str(np.sum(murs) / murs.size) + "\n"

    f.write(s)

    f.close()

    f = open('mur2x-eth.dat', 'a')

    s = str(w) + "  "
    for imur2 in range(murs2.size):
        s = s + str(murs2[imur2]) + "  "
    s = s + str(np.sum(murs2) / murs2.size) + "\n"

    f.write(s)

    f.close()

    iw += 1


text:
# coding: utf-8
# /*##########################################################################
# Copyright (C) 2016 European Synchrotron Radiation Facility
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
# ############################################################################*/

__authors__ = ["D. Naudet"]
__license__ = "MIT"
__date__ = "01/02/2016"

"""
Tests of the histogramnd function, error cases.
"""
import sys
import platform
import unittest

import numpy as np

from silx.math.chistogramnd import chistogramnd as histogramnd
from silx.math import Histogramnd

# ==============================================================
# ==============================================================
# ==============================================================

class _Test_chistogramnd_errors(unittest.TestCase):
    """
    Unit tests of the chistogramnd error cases.
    """
    __test__ = False  # ignore abstract class

    def setUp(self):
        self.skipTest("Abstract class")

    def test_weights_shape(self):
        """
        """

        for err_w_shape in self.err_weights_shapes:
            test_msg = ('Testing invalid weights shape : {0}'
                        ''.format(err_w_shape))

            err_weights = np.random.randint(0,
                                                    high=10,
                                                    size=err_w_shape)
            err_weights = err_weights.astype(np.double)

            ex_str = None
            try:
                histo, cumul = histogramnd(self.sample,
                                           self.histo_range,
                                           self.n_bins,
                                           weights=err_weights)[0:2]
            except ValueError as ex:
                ex_str = str(ex)

            self.assertIsNotNone(ex_str, msg=test_msg)
            self.assertEqual(ex_str,
                             '<weights> must be an array whose length '
                             'is equal to the number of samples.')

    def test_histo_range_shape(self):
        """
        """
        n_dims = 1 if len(self.s_shape) == 1 else self.s_shape[1]
        expected_txt_tpl = ('<histo_range> error : expected {n_dims} sets '
                            'of lower and upper bin edges, '
                            'got the following instead : {histo_range}. '
                            '(provided <sample> contains '
                            '{n_dims}D values)')

        for err_histo_range in self.err_histo_range_shapes:
            test_msg = ('Testing invalid histo_range shape : {0}'
                        ''.format(err_histo_range))

            expected_txt = expected_txt_tpl.format(histo_range=err_histo_range,
                                                   n_dims=n_dims)

            ex_str = None
            try:
                histo, cumul = histogramnd(self.sample,
                                           err_histo_range,
                                           self.n_bins,
                                           weights=self.weights)[0:2]
            except ValueError as ex:
                ex_str = str(ex)

            self.assertIsNotNone(ex_str, msg=test_msg)
            self.assertEqual(ex_str, expected_txt, msg=test_msg)

    def test_nbins_shape(self):
        """
        """

        expected_txt = ('n_bins must be either a scalar (same number '
                        'of bins for all dimensions) or '
                        'an array (number of bins for each '
                        'dimension).')

        for err_n_bins in self.err_n_bins_shapes:
            test_msg = ('Testing invalid n_bins shape : {0}'
                        ''.format(err_n_bins))

            ex_str = None
            try:
                histo, cumul = histogramnd(self.sample,
                                           self.histo_range,
                                           err_n_bins,
                                           weights=self.weights)[0:2]
            except ValueError as ex:
                ex_str = str(ex)

            self.assertIsNotNone(ex_str, msg=test_msg)
            self.assertEqual(ex_str, expected_txt, msg=test_msg)

    def test_nbins_values(self):
        """
        """
        expected_txt = ('<n_bins> : only positive values allowed.')

        for err_n_bins in self.err_n_bins_values:
            test_msg = ('Testing invalid n_bins value : {0}'
                        ''.format(err_n_bins))

            ex_str = None
            try:
                histo, cumul = histogramnd(self.sample,
                                           self.histo_range,
                                           err_n_bins,
                                           weights=self.weights)[0:2]
            except ValueError as ex:
                ex_str = str(ex)

            self.assertIsNotNone(ex_str, msg=test_msg)
            self.assertEqual(ex_str, expected_txt, msg=test_msg)

    def test_histo_shape(self):
        """
        """
        for err_h_shape in self.err_histo_shapes:

            # windows & python 2.7 : numpy shapes are long values
            if platform.system() == 'Windows':
                version = (sys.version_info.major, sys.version_info.minor)
                if version <= (2, 7):
                    err_h_shape = tuple([long(val) for val in err_h_shape])

            test_msg = ('Testing invalid histo shape : {0}'
                        ''.format(err_h_shape))

            expected_txt = ('Provided <histo> array doesn\'t have '
                            'a shape compatible with <n_bins> '
                            ': should be {0} instead of {1}.'
                            ''.format(self.h_shape, err_h_shape))

            histo = np.zeros(shape=err_h_shape, dtype=np.uint32)

            ex_str = None
            try:
                histo, cumul = histogramnd(self.sample,
                                           self.histo_range,
                                           self.n_bins,
                                           weights=self.weights,
                                           histo=histo)[0:2]
            except ValueError as ex:
                ex_str = str(ex)

            self.assertIsNotNone(ex_str, msg=test_msg)
            self.assertEqual(ex_str, expected_txt, msg=test_msg)

    def test_histo_dtype(self):
        """
        """
        for err_h_dtype in self.err_histo_dtypes:
            test_msg = ('Testing invalid histo dtype : {0}'
                        ''.format(err_h_dtype))

            histo = np.zeros(shape=self.h_shape, dtype=err_h_dtype)

            expected_txt = ('Provided <histo> array doesn\'t have '
                            'the expected type '
                            ': should be {0} instead of {1}.'
                            ''.format(np.uint32, histo.dtype))

            ex_str = None
            try:
                histo, cumul = histogramnd(self.sample,
                                           self.histo_range,
                                           self.n_bins,
                                           weights=self.weights,
                                           histo=histo)[0:2]
            except ValueError as ex:
                ex_str = str(ex)

            self.assertIsNotNone(ex_str, msg=test_msg)
            self.assertEqual(ex_str, expected_txt, msg=test_msg)

    def test_weighted_histo_shape(self):
        """
        """
        # using the same values as histo
        for err_h_shape in self.err_histo_shapes:

            # windows & python 2.7 : numpy shapes are long values
            if platform.system() == 'Windows':
                version = (sys.version_info.major, sys.version_info.minor)
                if version <= (2, 7):
                    err_h_shape = tuple([long(val) for val in err_h_shape])

            test_msg = ('Testing invalid weighted_histo shape : {0}'
                        ''.format(err_h_shape))

            expected_txt = ('Provided <weighted_histo> array doesn\'t have '
                            'a shape compatible with <n_bins> '
                            ': should be {0} instead of {1}.'
                            ''.format(self.h_shape, err_h_shape))

            cumul = np.zeros(shape=err_h_shape, dtype=np.double)

            ex_str = None
            try:
                histo, cumul = histogramnd(self.sample,
                                           self.histo_range,
                                           self.n_bins,
                                           weights=self.weights,
                                           weighted_histo=cumul)[0:2]
            except ValueError as ex:
                ex_str = str(ex)

            self.assertIsNotNone(ex_str, msg=test_msg)
            self.assertEqual(ex_str, expected_txt, msg=test_msg)

    def test_cumul_dtype(self):
        """
        """
        # using the same values as histo
        for err_h_dtype in self.err_histo_dtypes:
            test_msg = ('Testing invalid weighted_histo dtype : {0}'
                        ''.format(err_h_dtype))

            cumul = np.zeros(shape=self.h_shape, dtype=err_h_dtype)

            expected_txt = ('Provided <weighted_histo> array doesn\'t have '
                            'the expected type '
                            ': should be {0} or {1} instead of {2}.'
                            ''.format(np.float64, np.float32, cumul.dtype))

            ex_str = None
            try:
                histo, cumul = histogramnd(self.sample,
                                           self.histo_range,
                                           self.n_bins,
                                           weights=self.weights,
                                           weighted_histo=cumul)[0:2]
            except ValueError as ex:
                ex_str = str(ex)

            self.assertIsNotNone(ex_str, msg=test_msg)
            self.assertEqual(ex_str, expected_txt, msg=test_msg)

    def test_wh_histo_dtype(self):
        """
        """
        # using the same values as histo
        for err_h_dtype in self.err_histo_dtypes:
            test_msg = ('Testing invalid wh_dtype dtype : {0}'
                        ''.format(err_h_dtype))

            expected_txt = ('<wh_dtype> type not supported : {0}.'
                            ''.format(err_h_dtype))

            ex_str = None
            try:
                histo, cumul = histogramnd(self.sample,
                                           self.histo_range,
                                           self.n_bins,
                                           weights=self.weights,
                                           wh_dtype=err_h_dtype)[0:2]
            except ValueError as ex:
                ex_str = str(ex)

            self.assertIsNotNone(ex_str, msg=test_msg)
            self.assertEqual(ex_str, expected_txt, msg=test_msg)

    def test_unmanaged_dtypes(self):
        """
        """
        for err_unmanaged_dtype in self.err_unmanaged_dtypes:
            test_msg = ('Testing unmanaged dtypes : {0}'
                        ''.format(err_unmanaged_dtype))

            sample = self.sample.astype(err_unmanaged_dtype[0])
            weights = self.weights.astype(err_unmanaged_dtype[1])

            expected_txt = ('Case not supported - sample:{0} '
                            'and weights:{1}.'
                            ''.format(sample.dtype,
                                      weights.dtype))

            ex_str = None
            try:
                histogramnd(sample,
                            self.histo_range,
                            self.n_bins,
                            weights=weights)
            except TypeError as ex:
                ex_str = str(ex)

            self.assertIsNotNone(ex_str, msg=test_msg)
            self.assertEqual(ex_str, expected_txt, msg=test_msg)

    def test_uncontiguous_histo(self):
        """
        """
        # non contiguous array
        shape = np.array(self.n_bins, ndmin=1)
        shape[0] *= 2
        histo_tmp = np.zeros(shape)
        histo = histo_tmp[::2, ...]

        expected_txt = ('<histo> must be a C_CONTIGUOUS numpy array.')

        ex_str = None
        try:
            histogramnd(self.sample,
                        self.histo_range,
                        self.n_bins,
                        weights=self.weights,
                        histo=histo)
        except ValueError as ex:
            ex_str = str(ex)

        self.assertIsNotNone(ex_str)
        self.assertEqual(ex_str, expected_txt)

    def test_uncontiguous_weighted_histo(self):
        """
        """
        # non contiguous array
        shape = np.array(self.n_bins, ndmin=1)
        shape[0] *= 2
        cumul_tmp = np.zeros(shape)
        cumul = cumul_tmp[::2, ...]

        expected_txt = ('<weighted_histo> must be a C_CONTIGUOUS numpy array.')

        ex_str = None
        try:
            histogramnd(self.sample,
                        self.histo_range,
                        self.n_bins,
                        weights=self.weights,
                        weighted_histo=cumul)
        except ValueError as ex:
            ex_str = str(ex)

        self.assertIsNotNone(ex_str)
        self.assertEqual(ex_str, expected_txt)

class Test_chistogramnd_1D_errors(_Test_chistogramnd_errors):
    """
    Unit tests of the 1D histogramnd error cases.
    """
    __test__ = True  # because _Test_chistogramnd_errors is ignored

    def setUp(self):
        # nominal values
        self.n_elements = 1000
        self.s_shape = (self.n_elements,)
        self.w_shape = (self.n_elements,)

        self.histo_range = [0., 100.]
        self.n_bins = 10

        self.h_shape = (self.n_bins,)

        self.sample = np.random.randint(0,
                                                high=10,
                                                size=self.s_shape)
        self.sample = self.sample.astype(np.double)

        self.weights = np.random.randint(0,
                                                 high=10,
                                                 size=self.w_shape)
        self.weights = self.weights.astype(np.double)

        self.err_weights_shapes = ((self.n_elements+1,),
                                   (self.n_elements-1,),
                                   (self.n_elements-1, 3))
        self.err_histo_range_shapes = ([0.],
                                    [0., 1., 2.],
                                    [[0.], [1.]])
        self.err_n_bins_shapes = ([10, 2],
                                  [[10], [2]])
        self.err_n_bins_values = (0,
                                  [-10],
                                  None)
        self.err_histo_shapes = ((self.n_bins+1,),
                                 (self.n_bins-1,),
                                 (self.n_bins, self.n_bins))
        # these are used for testing the histo parameter as well
        #   as the weighted_histo parameter.
        self.err_histo_dtypes = (np.uint16,
                                 np.float16)

        self.err_unmanaged_dtypes = ((np.double, np.uint16),
                                     (np.uint16, np.double),
                                     (np.uint16, np.uint16))

class Test_chistogramnd_ND_range(unittest.TestCase):
    """

    """

    def test_invalid_histo_range(self):
        data = np.random.random((60, 60))
        nbins = 10

        with self.assertRaises(ValueError):
            histo_range = data.min(), np.inf

            Histogramnd(sample=data.ravel(),
                        histo_range=histo_range,
                        n_bins=nbins)

            histo_range = data.min(), np.nan

            Histogramnd(sample=data.ravel(),
                        histo_range=histo_range,
                        n_bins=nbins)

class Test_chistogramnd_ND_errors(_Test_chistogramnd_errors):
    """
    Unit tests of the 3D histogramnd error cases.
    """
    __test__ = True  # because _Test_chistogramnd_errors is ignored

    def setUp(self):
        # nominal values
        self.n_elements = 1000
        self.s_shape = (self.n_elements, 3)
        self.w_shape = (self.n_elements,)

        self.histo_range = [[0., 100.], [0., 100.], [0., 100.]]
        self.n_bins = (10, 20, 30)

        self.h_shape = self.n_bins

        self.sample = np.random.randint(0,
                                                high=10,
                                                size=self.s_shape)
        self.sample = self.sample.astype(np.double)

        self.weights = np.random.randint(0,
                                                 high=10,
                                                 size=self.w_shape)
        self.weights = self.weights.astype(np.double)

        self.err_weights_shapes = ((self.n_elements+1,),
                                   (self.n_elements-1,),
                                   (self.n_elements-1, 3))
        self.err_histo_range_shapes = ([0.],
                                    [0., 1.],
                                    [[0., 10.], [0., 10.]],
                                    [0., 10., 0, 10., 0, 10.])
        self.err_n_bins_shapes = ([10, 2],
                                  [[10], [20], [30]])
        self.err_n_bins_values = (0,
                                  [-10],
                                  [10, 20, -4],
                                  None,
                                  [10, None, 30])
        self.err_histo_shapes = ((self.n_bins[0]+1,
                                  self.n_bins[1],
                                  self.n_bins[2]),
                                 (self.n_bins[0],
                                  self.n_bins[1],
                                  self.n_bins[2]-1),
                                 (self.n_bins[0],
                                  self.n_bins[1]),
                                 (self.n_bins[1],
                                  self.n_bins[0],
                                  self.n_bins[2]),
                                 (self.n_bins[0],
                                  self.n_bins[1],
                                  self.n_bins[2],
                                  10)
                                 )
        # these are used for testing the histo parameter as well
        #   as the weighted_histo parameter.
        self.err_histo_dtypes = (np.uint16,
                                 np.float16)

        self.err_unmanaged_dtypes = ((np.double, np.uint16),
                                     (np.uint16, np.double),
                                     (np.uint16, np.uint16))

max_stars_repo_name: rnwatanabe/silx

max_issues_repo_name: rnwatanabe/silx

max_forks_repo_name: rnwatanabe/silx


text:
import numpy as np
from parallel_statistics import ParallelMeanVariance, ParallelMean

def read_shear_catalog_type(stage):
    """
    Determine the type of shear catalog a stage is using as input.
    Returns a string, e.g. metacal, lensfit.
    Also sets shear_catalog_type in the stage's configuration
    so that it is available later and is saved in output.
    """
    with stage.open_input('shear_catalog', wrapper=True) as f:
        shear_catalog_type = f.catalog_type
        stage.config['shear_catalog_type'] = shear_catalog_type
    return shear_catalog_type

def metacal_variants(*names):
    return [
        name + suffix
        for suffix in ['', '_1p', '_1m', '_2p', '_2m']
        for name in names
    ]

def band_variants(bands, *names, shear_catalog_type='metacal'):
    if shear_catalog_type=='metacal':
        return [
            name + "_" + band + suffix
            for suffix in ['', '_1p', '_1m', '_2p', '_2m']
            for band in bands
            for name in names
        ]
    else:
        return [
            name + "_" + band
            for band in bands
            for name in names
        ]

def calculate_selection_response(g1, g2, sel_1p, sel_2p, sel_1m, sel_2m, delta_gamma):
    import numpy as np

    S = np.ones((2,2))
    S_11 = (g1[sel_1p].mean() - g1[sel_1m].mean()) / delta_gamma
    S_12 = (g1[sel_2p].mean() - g1[sel_2m].mean()) / delta_gamma
    S_21 = (g2[sel_1p].mean() - g2[sel_1m].mean()) / delta_gamma
    S_22 = (g2[sel_2p].mean() - g2[sel_2m].mean()) / delta_gamma

    # Also save the selection biases as a matrix.
    S[0,0] = S_11
    S[0,1] = S_12
    S[1,0] = S_21
    S[1,1] = S_22

    return S

def calculate_shear_response(g1_1p,g1_2p,g1_1m,g1_2m,g2_1p,g2_2p,g2_1m,g2_2m,delta_gamma):
    import numpy as np

    n = len(g1_1p)
    R =  R = np.zeros((n,2,2))
    R_11 = (g1_1p - g1_1m) / delta_gamma
    R_12 = (g1_2p - g1_2m) / delta_gamma
    R_21 = (g2_1p - g2_1m) / delta_gamma
    R_22 = (g2_2p - g2_2m) / delta_gamma

    R[:,0,0] = R_11
    R[:,0,1] = R_12
    R[:,1,0] = R_21
    R[:,1,1] = R_22

    R = np.mean(R, axis=0)
    return R

def apply_metacal_response(R, S, g1, g2):
    # The values of R are assumed to already
    # have had appropriate weights included
    from numpy.linalg import pinv
    import numpy as np

    mcal_g = np.stack([g1,g2], axis=1)

    R_total = R+S

    # Invert the responsivity matrix
    Rinv = pinv(R_total)
    mcal_g = (Rinv @ mcal_g.T)

    return mcal_g[0], mcal_g[1]

def apply_lensfit_calibration(g1, g2, weight, c1=0, c2=0, sigma_e=0, m=0):
    w_tot = np.sum(weight)
    m = np.sum(weight*m)/w_tot        #if m not provided, default is m=0, so one_plus_K=1
    one_plus_K = 1.+m
    R = 1. - np.sum(weight*sigma_e)/w_tot
    g1 = (1./(one_plus_K))*((g1/R)-c1)
    g2 = (1./(one_plus_K))*((g2/R)-c2)
    return g1, g2, weight, one_plus_K

class _DataWrapper:
    """
    This little helper class wraps dictionaries
    or other mappings, and whenever something is
    looked up in it, it first checks if there is
    a column with the specified suffix present instead
    and returns that if so.
    """
    def __init__(self, data, suffix):
        """Create

        Parameters
        ----------
        data: dict

        suffix: str
        """
        self.suffix = suffix
        self.data = data

    def __getitem__(self, name):
        variant_name  = name + self.suffix
        if variant_name in self.data:
            return self.data[variant_name]
        else:
            return self.data[name]

    def __contains__(self, name):
        return (name in self.data)

class MetacalCalculator:
    """
    This class builds up the total response and selection calibration
    factors for Metacalibration from each chunk of data it is given.
    At the end an MPI communicator can be supplied to collect together
    the results from the different processes.

    To do this we need the function used to select the data, and the instance
    this function to each of the metacalibrated variants automatically by
    wrapping the data object passed in to it and modifying the names of columns
    that are looked up.
    """
    def __init__(self, selector, delta_gamma):
        """
        Initialize the Calibrator using the function you will use to select
        objects. That function should take at least one argument,
        the chunk of data to select on.  It should look up the original
        names of the columns to select on, without the metacal suffix.

        The MetacalCalculator will then wrap the data passed to it so that
        when a metacalibrated column is used for selection then the appropriate
        variant column is selected instead.

        The selector can take further *args and **kwargs, passed in when adding
        data.

        Parameters
        ----------
        selector: function
            Function that selects objects
        delta_gamma: float
            The difference in applied g between 1p and 1m metacal variants
        """
        self.selector = selector
        self.count = 0
        self.delta_gamma = delta_gamma
        self.cal_bias_means = ParallelMean(size=4)
        self.sel_bias_means = ParallelMean(size=8)

    def add_data(self, data, *args, **kwargs):
        """Select objects from a new chunk of data and tally their responses

        Parameters
        ----------
        data: dict
            Dictionary of data columns to select on and add
        *args
            Positional arguments to be passed to the selection function
        **kwargs
            Keyword arguments to be passed to the selection function

        """
        # These all wrap the catalog such that lookups find the variant
        # column if available.
        # For example, if I look up data_1p["x"] then it will check if
        # data["x_1p"] exists and return that if so.  Otherwise it will
        # fall back to "x"
        data_00 = _DataWrapper(data, '')
        data_1p = _DataWrapper(data, '_1p')
        data_1m = _DataWrapper(data, '_1m')
        data_2p = _DataWrapper(data, '_2p')
        data_2m = _DataWrapper(data, '_2m')

        # These are the selections from this chunk of data
        # that we would make under different shears, the baseline
        # and all the others.  self.selector is a function that the user
        # supplied in init, not a method
        sel_00 = self.selector(data_00, *args, **kwargs)
        sel_1p = self.selector(data_1p, *args, **kwargs)
        sel_1m = self.selector(data_1m, *args, **kwargs)
        sel_2p = self.selector(data_2p, *args, **kwargs)
        sel_2m = self.selector(data_2m, *args, **kwargs)

        g1 = data_00['mcal_g1']
        g2 = data_00['mcal_g2']
        weight = data_00['weight']
        n = g1[sel_00].size

        # Record the count for this chunk, for summation later
        self.count += n

        # This is the estimator response, correcting  bias of the shear estimator itself
        # We have four components, and want the weighted mean of each, which we use
        # the ParallelMean class to get
        w00 = weight[sel_00]
        R00 = data_1p['mcal_g1'][sel_00] - data_1m['mcal_g1'][sel_00]
        R01 = data_2p['mcal_g1'][sel_00] - data_2m['mcal_g1'][sel_00]
        R10 = data_1p['mcal_g2'][sel_00] - data_1m['mcal_g2'][sel_00]
        R11 = data_2p['mcal_g2'][sel_00] - data_2m['mcal_g2'][sel_00]

        # TODO: if there is a weight per variant would we use that here?
        # Not currently used though.
        self.cal_bias_means.add_data(0, R00, w00)
        self.cal_bias_means.add_data(1, R01, w00)
        self.cal_bias_means.add_data(2, R10, w00)
        self.cal_bias_means.add_data(3, R11, w00)

        # Now we handle the selection bias.  This value is given by the
        # difference in the (weighted) means of the shears under two selections,
        # the positive and negative shear, divided by the applied shear
        # See the second term in equation 14 of https://arxiv.org/pdf/1702.02601
        # We need to calculate eight means:
        # (g1 or g2) X (shear applied to 1 or 2) X (plus or minus shear)
        # We again use the ParallelMean class to handle this for us.

        self.sel_bias_means.add_data(0, g1[sel_1p], weight[sel_1p])
        self.sel_bias_means.add_data(1, g1[sel_1m], weight[sel_1m])
        self.sel_bias_means.add_data(2, g1[sel_2p], weight[sel_2p])
        self.sel_bias_means.add_data(3, g1[sel_2m], weight[sel_2m])
        self.sel_bias_means.add_data(4, g2[sel_1p], weight[sel_1p])
        self.sel_bias_means.add_data(5, g2[sel_1m], weight[sel_1m])
        self.sel_bias_means.add_data(6, g2[sel_2p], weight[sel_2p])
        self.sel_bias_means.add_data(7, g2[sel_2m], weight[sel_2m])

        # The user of this class may need the base selection, so return it
        return sel_00

    def collect(self, comm=None, allgather=False):
        """
        Finalize and sum up all the response values, returning separate
        R (estimator response) and S (selection bias) 2x2 matrices
        Parameters
        ----------
        comm: MPI Communicator
            If supplied, all processors response values will be combined together.
            All processes will return the same final value
        Returns
        -------
        R: 2x2 array
            Estimator response matrix
        S: 2x2 array
            Selection bias matrix
        """
        # collect all the things we need
        if comm is not None:
            if allgather:
                count = comm.allreduce(self.count)
            else:
                count = comm.reduce(self.count)
        else:
            count = self.count

        # Collect the mean values we need
        mode = ("allgather" if allgather else "gather")
        _, S = self.sel_bias_means.collect(comm, mode)
        _, R = self.cal_bias_means.collect(comm, mode)

        # Unpack the flat mean R and S values into
        # matrices
        R_mean = np.zeros((2, 2))
        R_mean[0, 0] = R[0]
        R_mean[0, 1] = R[1]
        R_mean[1, 0] = R[2]
        R_mean[1, 1] = R[3]
        R_mean /= self.delta_gamma

        # Need to take all the diffs to compute the
        # S term
        S_mean = np.zeros((2,2))
        S_mean[0, 0] = S[0] - S[1]
        S_mean[0, 1] = S[2] - S[3]
        S_mean[1, 0] = S[4] - S[5]
        S_mean[1, 1] = S[6] - S[7]
        S_mean /= self.delta_gamma

        return R_mean, S_mean, count

class LensfitCalculator:
    """
    This class builds up the total calibration
    factors for lensfit-convention shears from each chunk of data it is given.
    Note here we derive the c-terms from the data (in constrast to averaging
    values derived from simulations and stored in the catalog.)
    At the end an MPI communicator can be supplied to collect together
    the results from the different processes.

    """
    def __init__(self, selector, input_m_is_weighted=False):
        """
        Initialize the Calibrator using the function you will use to select
        objects. That function should take at least one argument,
        the chunk of data to select on.

        The selector can take further *args and **kwargs, passed in when adding
        data.

        Parameters
        ----------
        selector: function
            Function that selects objects
        """
        self.selector = selector
        # Create a set of calculators that will calculate (in parallel)
        # the three quantities we need to compute the overall calibration
        # We create these, then add data to them below, then collect them
        # together over all the processes
        self.K = ParallelMean(1)
        self.C = ParallelMean(2)
        self.count = 0
        self.input_m_is_weighted = input_m_is_weighted

    def add_data(self, data, *args, **kwargs):
        """Select objects from a new chunk of data and tally their responses

        Parameters
        ----------
        data: dict
            Dictionary of data columns to select on and add

        *args
            Positional arguments to be passed to the selection function
        **kwargs
            Keyword arguments to be passed to the selection function

        """
        # These all wrap the catalog such that lookups find the variant
        # column if available

        # This is just to let the selection tools access data.variant for feedback
        data = _DataWrapper(data, '')
        sel = self.selector(data, *args, **kwargs)

        # Extract the calibration quantities for the selected objects
        w = data['weight']
        K = data['m']
        g1 = data['g1']
        g2 = data['g2']
        n = g1[sel].size

         # Record the count for this chunk, for summation later
        self.count += n

        # Accumulate the calibration quantities so that later we
        # can compute the weighted mean of the values
        if self.input_m_is_weighted:
            # if the m values are already weighted don't use the weights here
            self.K.add_data(0, K[sel])
        else:
            # if not apply the weights
            self.K.add_data(0, K[sel], w[sel])
        self.C.add_data(0, g1[sel], w[sel])
        self.C.add_data(1, g2[sel], w[sel])

        return sel

    def collect(self, comm=None, allgather=False):
        """
        Finalize and sum up all the response values, returning calibration
        quantities.

        Parameters
        ----------
        comm: MPI Communicator
            If supplied, all processors response values will be combined together.
            All processes will return the same final value

        Returns
        -------

        K: float
            K = (1+m) calibration

        C: float array
            c1, c2 additive biases (weighted average of g1 and g2)

        """
        # The total number of objects is just the
        # number from all the processes summed together.
        if comm is not None:
            if allgather:
                count = comm.allreduce(self.count)
            else:
                count = comm.reduce(self.count)
        else:
            count = self.count

        # Collect the weighted means of these numbers.
        # this collects all the values from the different
        # processes and over all the chunks of data
        mode = ("allgather" if allgather else "gather")
        _ ,K = self.K.collect(comm, mode)
        _, C = self.C.collect(comm, mode)
        return K, C, count

class HSCCalculator:
    """
    This class builds up the total response calibration
    factors for HSC-convention shear-calibration from each chunk of data it is
    given.
    At the end an MPI communicator can be supplied to collect together
    the results from the different processes.

    """
    def __init__(self, selector):
        """
        Initialize the Calibrator using the function you will use to select
        objects. That function should take at least one argument,
        the chunk of data to select on.

        The selector can take further *args and **kwargs, passed in when adding
        data.

        Parameters
        ----------
        selector: function
            Function that selects objects
        """
        self.selector = selector
        # Create a set of calculators that will calculate (in parallel)
        # the three quantities we need to compute the overall calibration
        # We create these, then add data to them below, then collect them
        # together over all the processes
        self.K = ParallelMean(1)
        self.R = ParallelMean(1)
        self.count = 0

    def add_data(self, data, *args, **kwargs):
        """Select objects from a new chunk of data and tally their responses

        Parameters
        ----------
        data: dict
            Dictionary of data columns to select on and add

        *args
            Positional arguments to be passed to the selection function
        **kwargs
            Keyword arguments to be passed to the selection function

        """
        # These all wrap the catalog such that lookups find the variant
        # column if available

        # This is just to let the selection tools access data.variant for feedback
        data = _DataWrapper(data, '')
        sel = self.selector(data, *args, **kwargs)

        # Extract the calibration quantities for the selected objects
        w = data['weight']
        K = data['m']
        R = 1.0 - data['sigma_e'] ** 2
        n = w[sel].size
        self.count += w.size

        w = w[sel]

        # Accumulate the calibration quantities so that later we
        # can compute the weighted mean of the values
        self.R.add_data(0, R[sel], w)
        self.K.add_data(0, K[sel], w)

        return sel

    def collect(self, comm=None, allgather=False):
        """
        Finalize and sum up all the response values, returning calibration
        quantities.

        Parameters
        ----------
        comm: MPI Communicator
            If supplied, all processors response values will be combined together.
            All processes will return the same final value

        Returns
        -------
        R: float
            R calibration factor

        K: float
            K = (1+m) calibration

        N: int
            Total object count

        """
        # The total number of objects is just the
        # number from all the processes summed together.
        if comm is not None:
            if allgather:
                count = comm.allreduce(self.count)
            else:
                count = comm.reduce(self.count)
        else:
            count = self.count

        # Collect the weighted means of these numbers.
        # this collects all the values from the different
        # processes and over all the chunks of data
        mode = ("allgather" if allgather else "gather")
        _, R = self.R.collect(comm, mode)
        _ ,K = self.K.collect(comm, mode)
        return R, K, count

class MeanShearInBins:
    def __init__(self, x_name, limits, delta_gamma, cut_source_bin=False, shear_catalog_type='metacal'):
        self.x_name = x_name
        self.limits = limits
        self.delta_gamma = delta_gamma
        self.cut_source_bin = cut_source_bin
        self.shear_catalog_type = shear_catalog_type
        self.size = len(self.limits) - 1

        # We have to work out the mean g1, g2
        self.g1 = ParallelMeanVariance(self.size)
        self.g2 = ParallelMeanVariance(self.size)
        self.x  = ParallelMean(self.size)

        if shear_catalog_type=='metacal':
            self.calibrators = [MetacalCalculator(self.selector, delta_gamma) for i in range(self.size)]
        elif shear_catalog_type=='lensfit':
            self.calibrators = [LensfitCalculator(self.selector) for i in range(self.size)]
        elif shear_catalog_type=='hsc':
            self.calibrators = [HSCCalculator(self.selector) for i in range(self.size)]
        else:
            raise ValueError(f"Please specify metacal, lensfit or hsc for shear_catalog in config.")

    def selector(self, data, i):
        x = data[self.x_name]
        w = (x > self.limits[i]) & (x < self.limits[i+1])
        if self.cut_source_bin:
            w &= data['source_bin'] !=-1
        return np.where(w)

    def add_data(self, data):
        for i in range(self.size):
            w = self.calibrators[i].add_data(data, i)
            weight = data['weight'][w]
            if self.shear_catalog_type=='metacal':
                self.g1.add_data(i, data['mcal_g1'][w], weight)
                self.g2.add_data(i, data['mcal_g2'][w], weight)
            elif self.shear_catalog_type=='lensfit':
                self.g1.add_data(i, data['g1'][w], weight)
                self.g2.add_data(i, data['g2'][w], weight)
            elif self.shear_catalog_type=='hsc':
                self.g1.add_data(i, data['g1'][w]-data['c1'][w], weight)
                self.g2.add_data(i, data['g2'][w]-data['c2'][w], weight)
            self.x.add_data(i, data[self.x_name][w], weight)

    def collect(self, comm=None):
        count1, g1, var1 = self.g1.collect(comm, mode='gather')
        count2, g2, var2 = self.g2.collect(comm, mode='gather')

        _, mu = self.x.collect(comm, mode='gather')

        # Now we have the complete sample we can get the calibration matrix
        # to apply to it.
        R = []
        K = []
        C =[]
        for i in range(self.size):
            if self.shear_catalog_type=='metacal':
                # Tell the Calibrators to work out the responses
                r, s, _ = self.calibrators[i].collect(comm)
                # and record the total (a 2x2 matrix)
                R.append(r+s)
            elif self.shear_catalog_type=='lensfit':
                k, c, _ = self.calibrators[i].collect(comm)
                K.append(k)
                C.append(c)
            else:
                r, k, _ = self.calibrators[i].collect(comm)
                K.append(k)
                R.append(r)

        # Only the root processor does the rest
        if (comm is not None) and (comm.Get_rank() != 0):
            return None, None, None, None, None

        sigma1 = np.zeros(self.size)
        sigma2 = np.zeros(self.size)

        for i in range(self.size):
            # Get the shears and the errors on their means
            g = [g1[i], g2[i]]
            sigma = np.sqrt([var1[i]/count1[i], var2[i]/count2[i]])

            if self.shear_catalog_type=='metacal':
                # Get the inverse response matrix to apply
                R_inv = np.linalg.inv(R[i])

                # Apply the matrix in full to the shears and errors
                g1[i], g2[i] = R_inv @ g
                sigma1[i], sigma2[i] = R_inv @ sigma
            elif self.shear_catalog_type=='lensfit':
                g1[i] = g1[i]*(1./(1+K[i]))
                g2[i] = g2[i]*(1./(1+K[i]))

                sigma1[i] = (1./(1+K[i]))*(sigma[0])
                sigma2[i] = (1./(1+K[i]))*(sigma[1])
            else:
                g1[i] = (g1[i] / (2 * R[i]))/ (1 + K[i])
                g2[i] = (g2[i] / (2 * R[i]))/ (1 + K[i])

                sigma1[i] = (sigma[0] / (2 * R[i]))/ (1 + K[i])
                sigma2[i] = (sigma[1] / (2 * R[i]))/ (1 + K[i])

        return mu, g1, g2, sigma1, sigma2

max_stars_repo_name: Lhior/TXPipe

max_issues_repo_name: Lhior/TXPipe

max_forks_repo_name: Lhior/TXPipe

reason: import numpy,from numpy


text:
from typing import Any

import numpy as np

from .geocoder import Geocoder, ReverseGeocoder, NamesGeocoder

__all__ = [
    'distance',

    'geocode',
    'geocode_cities',
    'geocode_counties',
    'geocode_states',
    'geocode_countries',
    'reverse_geocode'
]

UNITS_DICT = {
    'mi': 3959,
    'km': 6371
}

GEOFUNC_TYPES = {
    'centroids': 'centroids',
    'boundaries': 'boundaries',
    'limits': 'limits',
    'region': 'regions'
}

def geocode(level=None, names=None, countries=None, states=None, counties=None, scope=None) -> NamesGeocoder:
    """
    Create a `NamesGeocoder`. Allows to refine ambiguous request with `where()` method,
    scope that limits area of geocoding or with parents.

    Parameters
    ----------
    level : {'country', 'state', 'county', 'city'}
        The level of administrative division. Autodetection by default.
    names : list or str
        Names of objects to be geocoded.
        For 'state' level: 'US-48' returns continental part of United States (48 states)
        in a compact form.
    countries : list
        Parent countries. Should have same size as names. Can contain strings or `Geocoder` objects.
    states : list
        Parent states. Should have same size as names. Can contain strings or `Geocoder` objects.
    counties : list
        Parent counties. Should have same size as names. Can contain strings or `Geocoder` objects.
    scope : str or `Geocoder`
        Limits area of geocoding. If parent country is set then error will be generated.
        If type is a string - geoobject should have geocoded scope in parents.
        If type is a `Geocoder` - geoobject should have geocoded scope in parents.
        Scope should contain only one entry.

    Returns
    -------
    `NamesGeocoder`
        Geocoder object specification.

    Examples
    --------
    .. jupyter-execute::
        :linenos:
        :emphasize-lines: 5

        from IPython.display import display
        from lets_plot import *
        from lets_plot.geo_data import *
        LetsPlot.setup_html()
        states = geocode('state').scope('Italy').get_boundaries(6)
        display(states.head())
        ggplot() + geom_map(data=states)

    |

    .. jupyter-execute::
        :linenos:
        :emphasize-lines: 5, 8

        from IPython.display import display
        from lets_plot import *
        from lets_plot.geo_data import *
        LetsPlot.setup_html()
        states = geocode(level='state', scope='US').get_geocodes()
        display(states.head())
        names = ['York'] * len(states.state)
        cities = geocode(names=names, states=states.state).ignore_not_found().get_centroids()
        display(cities.head())
        ggplot() + \\
            geom_livemap() + \\
            geom_point(data=cities, tooltips=layer_tooltips().line('@{found name}'))

    """
    return NamesGeocoder(level, names) \
        .scope(scope) \
        .countries(countries) \
        .states(states) \
        .counties(counties)

def geocode_cities(names=None) -> NamesGeocoder:
    """
    Create a `NamesGeocoder` object for cities. Allows to refine ambiguous request with
    `where()` method, with a scope that limits area of geocoding or with parents.

    Parameters
    ----------
    names : str or list
        Names of objects to be geocoded.

    Returns
    -------
    `NamesGeocoder`
        Geocoder object specification.

    Examples
    --------
    .. jupyter-execute::
        :linenos:
        :emphasize-lines: 5

        from IPython.display import display
        from lets_plot import *
        from lets_plot.geo_data import *
        LetsPlot.setup_html()
        cities = geocode_cities(['York', 'Jersey'])\\
                 .where(name='Jersey', scope='New Jersey').get_boundaries()
        display(cities)
        ggplot() + geom_map(aes(fill='found name'), data=cities, color='white')

    """
    return NamesGeocoder('city', names)

def geocode_counties(names=None) -> NamesGeocoder:
    """
    Create a `NamesGeocoder` object for counties. Allows to refine ambiguous request with
    `where()` method, with a scope that limits area of geocoding or with parents.

    Parameters
    ----------
    names : str or list
        Names of objects to be geocoded.

    Returns
    -------
    `NamesGeocoder`
        Geocoder object specification.

    Examples
    --------
    .. jupyter-execute::
        :linenos:
        :emphasize-lines: 5

        from IPython.display import display
        from lets_plot import *
        from lets_plot.geo_data import *
        LetsPlot.setup_html()
        counties = geocode_counties().scope('NY').get_boundaries(9)
        display(counties.head())
        ggplot() + geom_map(data=counties) + ggtitle('New York State Counties')

    """
    return NamesGeocoder('county', names)

def geocode_states(names=None) -> NamesGeocoder:
    """
    Create a `NamesGeocoder` object for states. Allows to refine ambiguous request with
    `where()` method, with a scope that limits area of geocoding or with parents.

    Parameters
    ----------
    names : str or list
        Names of objects to be geocoded.

    Returns
    -------
    `NamesGeocoder`
        Geocoder object specification.

    Examples
    --------
    .. jupyter-execute::
        :linenos:
        :emphasize-lines: 5

        from IPython.display import display
        from lets_plot import *
        from lets_plot.geo_data import *
        LetsPlot.setup_html()
        states = geocode_states().scope('UK').get_boundaries()
        display(states)
        ggplot() + \\
            geom_map(aes(fill='found name'), data=states, color='white') + \\
            ggtitle('UK States')

    """
    return NamesGeocoder('state', names)

def geocode_countries(names=None) -> NamesGeocoder:
    """
    Create a `NamesGeocoder` object for countries. Allows to refine ambiguous request with
    `where()` method.

    Parameters
    ----------
    names : str or list
        Names of objects to be geocoded.

    Returns
    -------
    `NamesGeocoder`
        Geocoder object specification.

    Examples
    --------
    .. jupyter-execute::
        :linenos:
        :emphasize-lines: 5

        from IPython.display import display
        from lets_plot import *
        from lets_plot.geo_data import *
        LetsPlot.setup_html()
        countries = geocode_countries(['Germany', 'Poland']).inc_res().get_boundaries()
        display(countries)
        ggplot() + geom_map(aes(fill='found name'), data=countries, color='white')

    """
    return NamesGeocoder('country', names)

def reverse_geocode(lon, lat, level=None, scope=None) -> ReverseGeocoder:
    """
    Convert a location as described by geographic coordinates to a `ReverseGeocoder` object.

    Parameters
    ----------
    lon : float
        Longitude coordinate of the geoobject.
    lat : float
        Latitude coordinate of the geoobject.
    level : {'country', 'state', 'county', 'city'}
        The level of administrative division.
    scope : str or `Geocoder`
        Specify this for resolving conflicts for disputed territories.

    Returns
    -------
    `ReverseGeocoder`
        Geocoder object specification.

    Examples
    --------
    .. jupyter-execute::
        :linenos:
        :emphasize-lines: 5

        from IPython.display import display
        from lets_plot import *
        from lets_plot.geo_data import *
        LetsPlot.setup_html()
        city = reverse_geocode(-73.87, 40.68, level='city').get_boundaries()
        display(city)
        ggplot() + geom_map(data=city) + ggtitle(city.iloc[0]['found name'])

    """
    return ReverseGeocoder(lon, lat, level, scope)

def distance(lon0, lat0, lon1, lat1, units='km'):
    """
    Calculate the distance between two points. Returns result in kilometers or miles.

    Parameters
    ----------
    lon0 : float
        Longitude coordinate of the first point.
    lat0 : float
        Latitude coordinate of the first point.
    lon1 : float
        Longitude coordinate of the second point.
    lat1 : float
        Latitude coordinate of the second point.
    units : {'mi', 'km'}, default='km'
        The units in which the result will be obtained.
        There are shorthands for values: 'mi' (miles), 'km' (kilometers).

    Returns
    -------
    float
        Distance between the points.

    Examples
    --------
    .. jupyter-execute::
        :linenos:
        :emphasize-lines: 3

        from lets_plot.geo_data import *
        cities = geocode_cities(['New York', 'Chicago']).get_centroids().geometry
        distance(cities[0].x, cities[0].y, cities[1].x, cities[1].y, units='mi')

    """
    return _calc_distance(lon0, lat0, lon1, lat1, units)

def _calc_distance(lon0, lat0, lon1, lat1, u):
    r = _prepare_units(u)

    lon0, lat0, lon1, lat1 = map(np.radians, [lon0, lat0, lon1, lat1])

    dlon = lon1 - lon0
    dlat = lat1 - lat0

    a = np.sin(dlat / 2.0) ** 2 + np.cos(lat0) * np.cos(lat1) * np.sin(dlon / 2.0) ** 2

    c = 2 * np.arcsin(np.sqrt(a))
    return c * r

def _prepare_units(units: Any) -> float:
    try:
        return UNITS_DICT[units]
    except KeyError:
        raise ValueError('Wrong units: {}. The units can take the following values: '
                         'mi (miles), km (kilometers).'.format(units))

max_stars_repo_name: IKrukov-HORIS/lets-plot

max_issues_repo_name: IKrukov-HORIS/lets-plot

max_forks_repo_name: IKrukov-HORIS/lets-plot


text:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
import seaborn as sns

def clean_df(df, col1, col2, look_for):
    '''
    INPUT:
    df - the pandas dataframe you want to search
    col1 - the column name you want to look through
    col2 - the column you want to count values from
    look_for - a list of strings you want to search for in each row of df[col]

    OUTPUT:
    new_df - a dataframe of each look_for with col2 values
    '''
    new_df = pd.DataFrame(columns=[col1,col2])
    #loop through list of ed types
    for val in look_for:
        #loop through rows
        for idx in range(df.shape[0]):
            #if the ed type is in the row add 1
            if df[col1][idx] == df[col1][idx] and df[col2][idx] == df[col2][idx]:
                if val in df[col1][idx]:
                    x2 = df[col2][idx]
                    x1 = val
                    new_df = new_df.append(pd.DataFrame([[x1,x2]],columns=[col1,col2]))
#     new_df = pd.DataFrame(pd.Series(new_df)).reset_index()
#     new_df.columns = [col1, col2]
#     new_df.sort_values('count', ascending=False, inplace=True)
    return new_df

##Function to group "At least once each week" and "At least once each day" together are 1

def frequency(row):
    if  row['StackOverflowFoundAnswer'] == "At least once each week" or row['StackOverflowFoundAnswer'] == "At least once each day":
        val = 1
    else:
        val = 0
    return val

def total_count(df, col1, col2, look_for):
    '''
    INPUT:
    df - the pandas dataframe you want to search
    col1 - the column name you want to look through
    col2 - the column you want to count values from
    look_for - a list of strings you want to search for in each row of df[col]

    OUTPUT:
    new_df - a dataframe of each look_for with the count of how often it shows up
    '''
    new_df = defaultdict(int)
    #loop through list of ed types
    for val in look_for:
        #loop through rows
        for idx in range(df.shape[0]):
            #if the ed type is in the row add 1
            if val in df[col1][idx]:
                new_df[val] += int(df[col2][idx])
    new_df = pd.DataFrame(pd.Series(new_df)).reset_index()
    new_df.columns = [col1, col2]
    new_df.sort_values('count', ascending=False, inplace=True)
    return new_df


text:
#################################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES), and is copyright (c) 2018-2021
# by the software owners: The Regents of the University of California, through
# Lawrence Berkeley National Laboratory,  National Technology & Engineering
# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia University
# Research Corporation, et al.  All rights reserved.
#
# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and
# license information.
#################################################################################
"""
Property data types.

Ability to import, etc. from text files is
part of the methods in the type.

Import property database from textfile(s):
* See :meth:`PropertyData.from_csv`, for the expected format for data.
* See :meth:`PropertyMetadata()` for the expected format for metadata.

"""
# stdlib
import csv
import json
import logging

# third-party
try:
    import pandas as pd
    import numpy as np
except ImportError:
    np, pd = None, None

# local
from .util import get_file
from . import tabular

__author__ = 'Dan Gunter'

_log = logging.getLogger(__name__)

class AddedCSVColumnError(KeyError):
    """Error for :meth:PropertyData.add_csv()
    """

    def __init__(self, names, how_bad, column_type=''):
        ctype = column_type + ' ' if column_type else ''
        if len(names) == 1:
            msg = 'Added CSV data {} {}column "{}"'.format(
                how_bad, ctype, list(names)[0]
            )
        else:
            msg = 'Added CSV data {} {}columns: {}'.format(
                how_bad, ctype, ', '.join(list(names))
            )
        KeyError.__init__(self, msg)

class Fields(tabular.Fields):
    """Constants for fields.
    """

    # Values for "type" field
    C_STATE, C_PROP = 'state', 'property'

class PropertyTable(tabular.Table):
    """Property data and metadata together (at last!)
    """

    def __init__(self, data=None, **kwargs):
        """Constructor.
        """
        if isinstance(data, PropertyData):
            pdata = data
        elif isinstance(data, list):
            pdata = PropertyData(data)
        else:
            raise TypeError('list or PropertyData object required')
        super(PropertyTable, self).__init__(data=pdata, **kwargs)

    @classmethod
    def load(cls, file_or_path, validate=True):
        """Create PropertyTable from JSON input.

        Args:
            file_or_path (file or str): Filename or file object
                from which to read the JSON-formatted data.
            validate (bool): If true, apply validation to input JSON data.

        Example input::

            {
                "meta": [
                    {"datatype": "MEA",
                     "info": "J. Chem. Eng. Data, 2009, Vol 54, pg. 306-310",
                     "notes": "r is MEA weight fraction in aqueous soln.",
                     "authors": "Amundsen, T.G., Lars, E.O., Eimer, D.A.",
                     "title": "Density and Viscosity of ..."}
                ],
                "data": [
                    {"name": "Viscosity Value",
                     "units": "mPa-s",
                     "values": [2.6, 6.2],
                     "error_type": "absolute",
                     "errors": [0.06, 0.004],
                     "type": "property"},
                    {"name": "r",
                     "units": "",
                     "values": [0.2, 1000],
                     "type": "state"}
                ]
            }
        """
        fp = get_file(file_or_path)
        d = json.load(fp)
        PropertyTable._validate_json(d)
        metalist = d[Fields.META]
        meta = [PropertyMetadata(m) for m in metalist]
        data = PropertyData(d[Fields.DATA])
        tbl = PropertyTable(data=data)
        for m in meta:
            tbl.add_metadata(m)
        return tbl

class PropertyData(tabular.TabularData):
    """Class representing property data that knows how to
    construct itself from a CSV file.

    You can build objects from multiple CSV files as well.
    See the property database section of the API docs for
    details, or read the code in :meth:`add_csv` and the
    tests in :mod:`idaes_dmf.propdb.tests.test_mergecsv`.
    """

    embedded_units = r'(.*)\((.*)\)'

    def __init__(self, data):
        """Construct new object from input list.

        Example input::

            [{
                "name": "Density Data",
                "units": "g/cm^3",
                "values": [1.0053, 1.0188, .., ],
                "errors": [.00005, .., .00005],
                "error_type": "absolute",
                "type": "property"
            }, ...etc...]

        Args:
            data (list): Input data columns

        Returns:
            (PropertyData) New instance.
        """
        super(PropertyData, self).__init__(data, error_column=True)
        self._nstates = len(self.states)

    @property
    def states(self):
        return [c for c in self.columns if self._is_state(c)]

    @property
    def properties(self):
        return [c for c in self.columns if self._is_prop(c)]

    @staticmethod
    def _is_state(c):
        return c[Fields.COLTYPE] == Fields.C_STATE

    @staticmethod
    def _is_prop(c):
        return c[Fields.COLTYPE] == Fields.C_PROP

    def names(self, states=True, properties=True):
        """Get column names.

        Args:
            states (bool): If False, exclude "state" data, e.g. the
                          ambient temperature, and only
                          include measured property values.
            properties (bool): If False, excluse property data

        Returns:
            list[str]: List of column names.
        """
        result = []
        if states:
            result.extend([v[Fields.DATA_NAME] for v in self.states])
        if properties:
            result.extend([v[Fields.DATA_NAME] for v in self.properties])
        return result

    def is_state_column(self, index):
        """Whether given column is state.

        Args:
            index (int): Index of column
        Returns:
            (bool) State or property and the column number.
        Raises:
            IndexError: No column at that index.
        """
        col = self.columns[index]
        return self._is_state(col)

    def is_property_column(self, index):
        """Whether given column is a property. See :meth:`is_state_column`."""
        return not self.is_state_column(index)

    def as_arr(self, states=True):
        """Export property data as arrays.

        Args:
            states (bool): If False, exclude "state" data, e.g. the
                          ambient temperature, and only
                          include measured property values.
        Returns:
            (values[M,N], errors[M,N]) Two arrays of floats,
            each with M columns having N values.
        Raises:
            ValueError if the columns are not all the same length
        """
        n, values, errors = None, [], []
        # extract state columns
        if states:
            for v in self.states:
                vals = v[Fields.DATA_VALUES]
                if n is None:
                    n = len(vals)
                elif len(vals) != n:
                    raise ValueError(
                        'State values "{}" length {} != {}'.format(
                            v[Fields.DATA_NAME], len(vals), n
                        )
                    )
                values.append(vals)
                errors.append([0] * len(vals))
        # extract property columns
        for v in self.properties:
            vals = v[Fields.DATA_VALUES]
            if n is None:
                n = len(vals)
            elif len(vals) != n:
                raise ValueError(
                    'Property values "{}" length {} != {}'.format(
                        v[Fields.DATA_NAME], len(vals), n
                    )
                )
            values.append(v[Fields.DATA_VALUES])
            errors.append(v[Fields.DATA_ERRORS])
        return values, errors

    def values_dataframe(self, states=True):
        """Get values as a dataframe.

        Args:
            states (bool): see :meth:`names()`.

        Returns:
            (pd.DataFrame) Pandas dataframe for values.

        Raises:
            ImportError: If `pandas` or `numpy` were never
                successfully imported.
        """
        return self._get_prop_dataframe(Fields.DATA_VALUES, states)

    def errors_dataframe(self, states=False):
        """Get errors as a dataframe.

        Args:
            states (bool): If False, exclude state data.
                           This is the default, because states do not
                           normally have associated error information.
        Returns:
            pd.DataFrame: Pandas dataframe for values.

        Raises:
            ImportError: If `pandas` or `numpy` were never
                successfully imported.
        """
        return self._get_prop_dataframe(Fields.DATA_ERRORS, states)

    def _get_prop_dataframe(self, field, states):
        self._check_pandas_import()
        a1, names = [], []
        if states:
            a1 = [v[field] for v in self.states]
            names = [v[Fields.DATA_NAME] for v in self.states]
        a1.extend([v[field] for v in self.properties])
        names.extend([v[Fields.DATA_NAME] for v in self.properties])
        a2 = np.array(a1).transpose()
        return pd.DataFrame(a2, columns=names)

    @staticmethod
    def from_csv(file_or_path, nstates=0):
        """Import the CSV data.

        Expected format of the  files is a header plus data rows.

        Header: Index-column,  Column-name(1), Error-column(1),  \
                Column-name(2), Error-column(2), ..
        Data: <index>, <val>, <errval>, <val>, <errval>, ..

        Column-name is in the format "Name (units)"

        Error-column is in the format "<type> Error", where "<type>" is
        the error type.

        Args:
            file_or_path (file-like or str): Input file
            nstates (int): Number of state columns, appearing
               first before property columns.

        Returns:
            PropertyData: New properties instance
        """
        input_file = get_file(file_or_path)
        csv_file = csv.reader(input_file)
        row = next(csv_file)
        names, data = PropertyData._prop_parse_csv_headers(nstates, row)
        for row in csv_file:
            # print('@@ parse csv row: {}'.format(row))
            PropertyData._parse_csv_row(data, row, error_column=True)
        obj = PropertyData(data)
        return obj

    def add_csv(self, file_or_path, strict=False):
        """Add to existing object from a new CSV file.

        Depending on the value of the `strict` argument (see
        below), the new file may or may not have the same
        properties as the object -- but it always needs to have
        the same number of state columns, and in the same order.

        .. note:: Data that is "missing" because of property columns in
           one CSV and not the other will be filled with `float(nan)` values.

        Args:
            file_or_path (file or str): Input file. This should be in exactly
                the same format as expected by :meth:from_csv().
            strict (bool): If true, require that the columns in the input
                CSV match columns in this object. Otherwise, only require
                that *state* columns in input CSV match columns in this
                object. New property columns are added, and matches
                to existing property columns will append the data.

        Raises:
            AddedCSVColumnError: If the new CSV column headers are not the
                same as the ones in this object.

        Returns:
            (int) Number of added rows
        """
        nstates = self._nstates
        input_file = get_file(file_or_path)
        csv_file = csv.reader(input_file)

        # Parse the header
        row = next(csv_file)
        hdr_names, hdr_data = PropertyData._prop_parse_csv_headers(nstates, row)

        # print('@@ add_csv, column names = {}, data columns = {}'
        #      .format(hdr_names, self.names()))

        # Check that set of keys in new data is the same
        cur_keys = set(self.names())
        new_keys = set(hdr_names)
        # This is used to re-order input data
        rowmap = None
        if strict:
            if cur_keys > new_keys:
                missing = cur_keys - new_keys
                raise AddedCSVColumnError(missing, 'is missing')
            elif new_keys > cur_keys:
                extra = new_keys - cur_keys
                raise AddedCSVColumnError(extra, 'has extra')
            elif new_keys != cur_keys:
                extra = new_keys - cur_keys
                missing = cur_keys - new_keys
                namelist = (
                    '(' + ','.join(extra) + ')',
                    'instead of',
                    '(' + ','.join(missing) + ')',
                )
                raise AddedCSVColumnError(namelist, 'has different')
        else:
            # check that all states are in common
            hdr_states = filter(self._is_state, hdr_data)
            new_states = [s[Fields.DATA_NAME] for s in hdr_states]
            new_states = set(new_states)
            cur_states = set(self.names(properties=False))
            if new_states != cur_states:
                extra = new_states - cur_states
                missing = cur_states - new_states
                if extra and missing:
                    namelist = (
                        '(' + ','.join(extra) + ')',
                        'instead of',
                        '(' + ','.join(missing) + ')',
                    )
                    raise AddedCSVColumnError(
                        namelist, 'has different', column_type='state'
                    )
                elif extra:
                    raise AddedCSVColumnError(extra, 'has extra', column_type='state')
                elif missing:
                    raise AddedCSVColumnError(
                        missing, 'is missing', column_type='state'
                    )
                else:
                    raise RuntimeError('unexpected branch')
            # check that at least one property is in common
            new_prop = new_keys - new_states
            if not new_prop:
                return 0  # no data
            cur_prop = set(self.names(states=False))

            # Add columns for all properties only found on the input,
            # and initialize values to a list of NaN's as long as the
            # current table, so data in all fields will be the same length.

            # Initialize rowmap with mapping for state columns
            rowmap = [-1] * len(hdr_names)
            idx = 0
            for i, s in enumerate(hdr_data):
                if s[Fields.COLTYPE] == Fields.C_PROP:
                    continue
                rowmap[i] = idx
                idx += 1

            nan_list = [float('nan')] * self.num_rows
            idx = 0
            for i, value in enumerate(hdr_data):
                if value[Fields.COLTYPE] == Fields.C_STATE:
                    continue
                name = value[Fields.DATA_NAME]
                if name not in cur_prop:
                    value[Fields.DATA_NAME] = name
                    value[Fields.DATA_VALUES] = nan_list[:]
                    value[Fields.DATA_ERRORS] = nan_list[:]
                    value[Fields.COLTYPE] = Fields.C_PROP
                    self._data.append(value)
                    rowmap[i] = len(self.properties) - 1
                else:
                    rowmap[i] = idx + self._nstates
                idx += 1

        # print("@@ rowmap = {}".format(rowmap))

        # Parse the new data
        num_added = 0
        new_rowlen = 1 + 2 * len(self.names())
        for row in csv_file:
            if rowmap:
                # Re-order according to the rowmap.
                # By initializing with NaN, any columns not in the
                # input, but in the current data, will be replaced with NaN
                # values.
                row2 = [float('nan')] * new_rowlen
                # print('@@ row={} row2-init={}'.format(row, row2))
                for i, j in enumerate(rowmap):
                    row2[j * 2 + 1] = row[i * 2 + 1]  # value
                    row2[j * 2 + 2] = row[i * 2 + 2]  # error
                row = row2
            self._parse_csv_row(self._data, row, error_column=True)
            num_added += 1
            self._nrows += 1

        return num_added

    @classmethod
    def _prop_parse_csv_headers(cls, nstates, headers):
        """Parse a row of CSV headers which are pairs
        of columns like "<name> [(units)], <error-type> Error".

        Returns:
             (names, data). Names is a list of all the column names.
             Data is a dict with two keys, "properties" and "states".
             Each value will be a list of property/state objects.
        """
        names, data = cls._parse_csv_headers(headers, error_column=True)
        for i in range(0, nstates):
            data[i][Fields.COLTYPE] = Fields.C_STATE
        for i in range(nstates, len(data)):
            data[i][Fields.COLTYPE] = Fields.C_PROP
        return names, data

class PropertyMetadata(tabular.Metadata):
    """Class to import property metadata.
    """

    pass

class PropertyColumn(tabular.Column):
    """Data column for a property.
    """

    type_name = 'Property'

    def __init__(self, name, data):
        tabular.Column.__init__(self, name, data)
        self.errors = data[Fields.DATA_ERRORS]
        self.error_type = data[Fields.DATA_ERRTYPE]

    def data(self):
        return {
            Fields.DATA_UNITS: self.units,
            Fields.DATA_VALUES: self.values,
            Fields.DATA_ERRORS: self.errors,
            Fields.DATA_ERRTYPE: self.error_type,
        }

class StateColumn(tabular.Column):
    """Data column for a state.
    """

    type_name = 'State'

    def __init__(self, name, data):
        tabular.Column.__init__(self, name, data)
        self.errors = [0.0] * len(self)
        self.error_type = 'none'

    def data(self):
        return {Fields.DATA_UNITS: self.units, Fields.DATA_VALUES: self.values}

def convert_csv(meta_csv, datatype, data_csv, nstates, output):
    meta = PropertyMetadata.from_csv(meta_csv)
    meta.datatype = datatype
    data = PropertyData.from_csv(data_csv, nstates)
    obj = PropertyTable(data=data, metadata=meta)
    ofile = get_file(output, mode='w')
    obj.dump(ofile)

max_stars_repo_name: carldlaird/idaes-pse

max_issues_repo_name: carldlaird/idaes-pse

max_forks_repo_name: carldlaird/idaes-pse


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

# In[ ]:

####################################<<<<Breast_cancer_prediction>>>>>>####################################

# In[ ]:

#part(1)--By:Manar Moeanse

# In[1]:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# In[2]:

DB = pd.read_csv('Breast_cancer_data.csv')
DB

# In[3]:

DB.head(5)

# In[4]:

DB.describe()

# In[ ]:

DB.info()

# In[ ]:

#part(2)--By:Mariam Mamdoh

# In[5]:

uneff = DB[DB.diagnosis == 0]
eff = DB[DB.diagnosis == 1]
len(uneff)

# In[ ]:

len(eff)

# In[6]:

uneffected = (len(uneff)/len(DB)) *100
print('people are uneffected = ', uneffected , '% .')
effected = (len(eff)/len(DB)) *100
print('people are effected = ', effected , '% .')

# In[ ]:

#part(3)--By:Hemat Shawky.

# In[7]:

plt.scatter(DB['diagnosis'],DB['mean_area'])

# In[8]:

plt.scatter(DB['mean_area'],DB['mean_texture'])

# In[9]:

plt.scatter(DB['mean_radius'],DB['mean_perimeter'])

# In[10]:

import seaborn as sns
sns.pairplot(data=DB)

# In[ ]:

#part(4)--By:Hager Mohamed.

# In[11]:

x = DB.drop('diagnosis', 1)
y = DB['diagnosis']
x

# In[12]:

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)
print(x_train.shape,x_test.shape)
print(y_train.shape,y_test.shape)

# In[13]:

from sklearn.linear_model import LinearRegression
model = LinearRegression ()

# In[14]:

model.fit(x_train,y_train)

# In[15]:

pred =model.predict(x_test)

# In[16]:

from sklearn.metrics import mean_squared_error

# In[17]:

error=np.sqrt(mean_squared_error(y_pred=pred,y_true=y_test))
print(error)

# In[18]:

print(model.score(x_test,y_test))

# In[ ]:

"""" BY:
        1-Manar Moeanse.
        2-Mariam Mamdoh.
        3-Hemat Shawky.
        4-Hager Mohamed.


text:
'''
### Recurrent Neural Network representation of C. elegans neural timeseries data
Script to model C. elegans neural timeseries data,
and generate similar data, using a Recurrent
Neural Network (specifically, a long short term memory or LSTM network)

Inspired by
http://karpathy.github.io/2015/05/21/rnn-effectiveness/

Based on code from
https://github.com/fchollet/keras/blob/master/examples/lstm_text_generation.py
'''

from keras.models import Sequential
from keras.layers import Layer
from keras.layers.core import Dense, Activation, Dropout
from keras.layers.recurrent import LSTM
from keras.utils.data_utils import get_file
from keras.optimizers import RMSprop

import numpy as np
import random
import sys
import os

import theano
from theano import tensor as T, function, printing

import biodatamanager as dm

os.environ["THEANO_FLAGS"] = "mode=FAST_RUN,device=gpu,floatX=float32"

# Layer implementing a Gaussian mixture model
# Implementation follows
# https://github.com/jrieke/lstm-biology/blob/master/prediction_stateful_mixture.ipynb
# and
# https://github.com/fchollet/keras/issues/1061
# Thanks to @jrieke for providing code and guidance

class GMMActivation(Layer):
    """
    GMM-like activation function.
    Assumes that input has (D+2)*M dimensions, where D is the dimensionality of the
    target data. The first M*D features are treated as means, the next M features as
    standard devs and the last M features as mixture components of the GMM.
    """
    def __init__(self, M, **kwargs):
        super(GMMActivation, self).__init__(**kwargs)
        self.M = M

    def get_output(self, train=False):
        X = self.get_input(train)
        D = T.shape(X)[1]/self.M - 2
        # leave mu values as they are since they're unconstrained
        # scale sigmas with exp, s.t. all values are non-negative
        X = T.set_subtensor(X[:,D*self.M:(D+1)*self.M], T.exp(X[:,D*self.M:(D+1)*self.M]))
        # scale alphas with softmax, s.t. that all values are between [0,1] and sum up to 1
        X = T.set_subtensor(X[:,(D+1)*self.M:(D+2)*self.M], T.nnet.softmax(X[:,(D+1)*self.M:(D+2)*self.M]))
        return X

    def get_config(self):
        config = {"name": self.__class__.__name__,
                  "M": self.M}
        base_config = super(GMMActivation, self).get_config()
        return dict(list(base_config.items()) + list(config.items()))

def gmm_loss(y_true, y_pred):
    """
    GMM loss function.
    Assumes that y_pred has (D+2)*M dimensions and y_true has D dimensions. The first
    M*D features are treated as means, the next M features as standard devs and the last
    M features as mixture components of the GMM.
    """
    def loss(m, M, D, y_true, y_pred):
        mu = y_pred[:,D*m:(m+1)*D]
        sigma = y_pred[:,D*M+m]
        alpha = y_pred[:,(D+1)*M+m]

        return (alpha/sigma/np.sqrt(2. * np.pi)) * T.exp(-T.sum(T.sqr(mu-y_true),-1)/(2*sigma**2))

    D = T.shape(y_true)[1]
    M = T.shape(y_pred)[1]/(D+2)
    seq = T.arange(M)
    result, _ = theano.scan(fn=loss, outputs_info=None,
    sequences=seq, non_sequences=[M, D, y_true, y_pred])
    return -T.log(result.sum(0))

def generate(seq, maxlen=1, bs=3, ep=5, output_iterations=10, num_mixture_components=10):
    # seq is a single sample, in the format (timesteps, features) !
    # TODO: expand code to support multiple samples, fed into model together as a batch
    # Cut the timeseries data (variable name 'seq') into semi-redundant sequence chunks of maxlen

    X = []
    y = []

    for i in range(0, len(seq) - maxlen):
        X.append(seq[i:i+maxlen])
        y.append(seq[i+maxlen])

    dim = len((X[0][0]))

    print("sequence chunks:", len(X))
    print("chunk width:", len(X[0]))
    print("vector dimension:", dim)
    print("number of mixture components:", num_mixture_components)

    X = np.array(X)
    y = np.array(y)

    # build the model: 2 stacked LSTM
    print('Build model...')
    model = Sequential()
    model.reset_states()
    model.add(LSTM((dim+2) * num_mixture_components, return_sequences=False, input_shape=(maxlen, dim)))
    model.add(Dense((dim+2) * num_mixture_components))

    model.add(GMMActivation(num_mixture_components))

    model.compile(loss=gmm_loss, optimizer=RMSprop(lr=0.001))

    # Train the model
    model.fit(X, y, batch_size=bs, nb_epoch=ep)

    # Generate timeseries
    x_seed = X[len(X)-1] #choose final in-sample data point to initialize model
    x_array = []
    x_array.append(x_seed)
    x = np.array(x_array)

    predicted = []
    for i in range(output_iterations):
        pred_parameters = model.predict_on_batch(x)[0]

        means = pred_parameters[:num_mixture_components * dim]
        sds = pred_parameters[(num_mixture_components * dim):(num_mixture_components * (dim+1))]
        weights = pred_parameters[(num_mixture_components * (dim + 1)):]

        print(means)
        print(sds)
        print(weights)

        means = means.reshape(num_mixture_components, dim)
        sds = sds[:, np.newaxis]
        weights = weights[:, np.newaxis]

        pred = weights * np.random.normal(means, sds)
        pred = np.sum(pred, axis=0)
        predicted.append(pred)

    return predicted

max_stars_repo_name: openworm/neuronal-analysis

max_issues_repo_name: openworm/neuronal-analysis

max_forks_repo_name: openworm/neuronal-analysis

reason: import numpy,import theano,from theano