text:
"""MXNet modules for graph convolutions."""
# pylint: disable= no-member, arguments-differ, invalid-name
import math
import mxnet as mx
from mxnet import gluon, nd
from mxnet.gluon import nn
import numpy as np

from . import utils
from ... import function as fn

__all__ = ['GraphConv', 'TAGConv', 'RelGraphConv']

class GraphConv(gluon.Block):
    r"""Apply graph convolution over an input signal.

    Graph convolution is introduced in `GCN <https://arxiv.org/abs/1609.02907>`__
    and can be described as below:

    .. math::
      h_i^{(l+1)} = \sigma(b^{(l)} + \sum_{j\in\mathcal{N}(i)}\frac{1}{c_{ij}}h_j^{(l)}W^{(l)})

    where :math:`\mathcal{N}(i)` is the neighbor set of node :math:`i`. :math:`c_{ij}` is equal
    to the product of the square root of node degrees:
    :math:`\sqrt{|\mathcal{N}(i)|}\sqrt{|\mathcal{N}(j)|}`. :math:`\sigma` is an activation
    function.

    The model parameters are initialized as in the
    `original implementation <https://github.com/tkipf/gcn/blob/master/gcn/layers.py>`__ where
    the weight :math:`W^{(l)}` is initialized using Glorot uniform initialization
    and the bias is initialized to be zero.

    Notes
    -----
    Zero in degree nodes could lead to invalid normalizer. A common practice
    to avoid this is to add a self-loop for each node in the graph, which
    can be achieved by:

    >>> g = ... # some DGLGraph
    >>> g.add_edges(g.nodes(), g.nodes())

    Parameters
    ----------
    in_feats : int
        Number of input features.
    out_feats : int
        Number of output features.
    norm : bool, optional
        If True, the normalizer :math:`c_{ij}` is applied. Default: ``True``.
    bias : bool, optional
        If True, adds a learnable bias to the output. Default: ``True``.
    activation: callable activation function/layer or None, optional
        If not None, applies an activation function to the updated node features.
        Default: ``None``.

    Attributes
    ----------
    weight : mxnet.gluon.parameter.Parameter
        The learnable weight tensor.
    bias : mxnet.gluon.parameter.Parameter
        The learnable bias tensor.
    """
    def __init__(self,
                 in_feats,
                 out_feats,
                 norm=True,
                 bias=True,
                 activation=None):
        super(GraphConv, self).__init__()
        self._in_feats = in_feats
        self._out_feats = out_feats
        self._norm = norm

        with self.name_scope():
            self.weight = self.params.get('weight', shape=(in_feats, out_feats),
                                          init=mx.init.Xavier(magnitude=math.sqrt(2.0)))
            if bias:
                self.bias = self.params.get('bias', shape=(out_feats,),
                                            init=mx.init.Zero())
            else:
                self.bias = None

        self._activation = activation

    def forward(self, graph, feat):
        r"""Compute graph convolution.

        Notes
        -----
            * Input shape: :math:`(N, *, \text{in_feats})` where * means any number of additional
              dimensions, :math:`N` is the number of nodes.
            * Output shape: :math:`(N, *, \text{out_feats})` where all but the last dimension are
              the same shape as the input.

        Parameters
        ----------
        graph : DGLGraph
            The graph.
        feat : mxnet.NDArray
            The input feature

        Returns
        -------
        mxnet.NDArray
            The output feature
        """
        graph = graph.local_var()
        if self._norm:
            degs = graph.in_degrees().astype('float32')
            norm = mx.nd.power(mx.nd.clip(degs, a_min=1, a_max=float("inf")), -0.5)
            shp = norm.shape + (1,) * (feat.ndim - 1)
            norm = norm.reshape(shp).as_in_context(feat.context)
            feat = feat * norm

        if self._in_feats > self._out_feats:
            # mult W first to reduce the feature size for aggregation.
            feat = mx.nd.dot(feat, self.weight.data(feat.context))
            graph.ndata['h'] = feat
            graph.update_all(fn.copy_src(src='h', out='m'),
                             fn.sum(msg='m', out='h'))
            rst = graph.ndata.pop('h')
        else:
            # aggregate first then mult W
            graph.ndata['h'] = feat
            graph.update_all(fn.copy_src(src='h', out='m'),
                             fn.sum(msg='m', out='h'))
            rst = graph.ndata.pop('h')
            rst = mx.nd.dot(rst, self.weight.data(feat.context))

        if self._norm:
            rst = rst * norm

        if self.bias is not None:
            rst = rst + self.bias.data(rst.context)

        if self._activation is not None:
            rst = self._activation(rst)

        return rst

    def __repr__(self):
        summary = 'GraphConv('
        summary += 'in={:d}, out={:d}, normalization={}, activation={}'.format(
            self._in_feats, self._out_feats,
            self._norm, self._activation)
        summary += '\n)'
        return summary

class TAGConv(gluon.Block):
    r"""Apply Topology Adaptive Graph Convolutional Network

    .. math::
        \mathbf{X}^{\prime} = \sum_{k=0}^K \mathbf{D}^{-1/2} \mathbf{A}
        \mathbf{D}^{-1/2}\mathbf{X} \mathbf{\Theta}_{k},

    where :math:`\mathbf{A}` denotes the adjacency matrix and
    :math:`D_{ii} = \sum_{j=0} A_{ij}` its diagonal degree matrix.

    Parameters
    ----------
    in_feats : int
        Number of input features.
    out_feats : int
        Number of output features.
    k: int, optional
        Number of hops :math: `k`. (default: 2)
    bias: bool, optional
        If True, adds a learnable bias to the output. Default: ``True``.
    activation: callable activation function/layer or None, optional
        If not None, applies an activation function to the updated node features.
        Default: ``None``.

    Attributes
    ----------
    lin : mxnet.gluon.parameter.Parameter
        The learnable weight tensor.
    bias : mxnet.gluon.parameter.Parameter
        The learnable bias tensor.
    """
    def __init__(self,
                 in_feats,
                 out_feats,
                 k=2,
                 bias=True,
                 activation=None):
        super(TAGConv, self).__init__()
        self.out_feats = out_feats
        self.k = k
        self.bias = bias
        self.activation = activation
        self.in_feats = in_feats

        self.lin = self.params.get(
            'weight', shape=(self.in_feats * (self.k + 1), self.out_feats),
            init=mx.init.Xavier(magnitude=math.sqrt(2.0)))
        if self.bias:
            self.h_bias = self.params.get('bias', shape=(out_feats,),
                                          init=mx.init.Zero())

    def forward(self, graph, feat):
        r"""Compute graph convolution

        Parameters
        ----------
        graph : DGLGraph
            The graph.
        feat : mxnet.NDArray
            The input feature of shape :math:`(N, D_{in})` where :math:`D_{in}`
            is size of input feature, :math:`N` is the number of nodes.

        Returns
        -------
        mxnet.NDArray
            The output feature of shape :math:`(N, D_{out})` where :math:`D_{out}`
            is size of output feature.
        """
        graph = graph.local_var()

        degs = graph.in_degrees().astype('float32')
        norm = mx.nd.power(mx.nd.clip(degs, a_min=1, a_max=float("inf")), -0.5)
        shp = norm.shape + (1,) * (feat.ndim - 1)
        norm = norm.reshape(shp).as_in_context(feat.context)

        rst = feat
        for _ in range(self.k):
            rst = rst * norm
            graph.ndata['h'] = rst

            graph.update_all(fn.copy_src(src='h', out='m'),
                             fn.sum(msg='m', out='h'))
            rst = graph.ndata['h']
            rst = rst * norm
            feat = mx.nd.concat(feat, rst, dim=-1)

        rst = mx.nd.dot(feat, self.lin.data(feat.context))
        if self.bias is not None:
            rst = rst + self.h_bias.data(rst.context)

        if self.activation is not None:
            rst = self.activation(rst)

        return rst

class RelGraphConv(gluon.Block):
    r"""Relational graph convolution layer.

    Relational graph convolution is introduced in "`Modeling Relational Data with Graph
    Convolutional Networks <https://arxiv.org/abs/1703.06103>`__"
    and can be described as below:

    .. math::

      h_i^{(l+1)} = \sigma(\sum_{r\in\mathcal{R}}
      \sum_{j\in\mathcal{N}^r(i)}\frac{1}{c_{i,r}}W_r^{(l)}h_j^{(l)}+W_0^{(l)}h_i^{(l)})

    where :math:`\mathcal{N}^r(i)` is the neighbor set of node :math:`i` w.r.t. relation
    :math:`r`. :math:`c_{i,r}` is the normalizer equal
    to :math:`|\mathcal{N}^r(i)|`. :math:`\sigma` is an activation function. :math:`W_0`
    is the self-loop weight.

    The basis regularization decomposes :math:`W_r` by:

    .. math::

      W_r^{(l)} = \sum_{b=1}^B a_{rb}^{(l)}V_b^{(l)}

    where :math:`B` is the number of bases.

    The block-diagonal-decomposition regularization decomposes :math:`W_r` into :math:`B`
    number of block diagonal matrices. We refer :math:`B` as the number of bases.

    Parameters
    ----------
    in_feat : int
        Input feature size.
    out_feat : int
        Output feature size.
    num_rels : int
        Number of relations.
    regularizer : str
        Which weight regularizer to use "basis" or "bdd"
    num_bases : int, optional
        Number of bases. If is none, use number of relations. Default: None.
    bias : bool, optional
        True if bias is added. Default: True
    activation : callable, optional
        Activation function. Default: None
    self_loop : bool, optional
        True to include self loop message. Default: False
    dropout : float, optional
        Dropout rate. Default: 0.0
    """
    def __init__(self,
                 in_feat,
                 out_feat,
                 num_rels,
                 regularizer="basis",
                 num_bases=None,
                 bias=True,
                 activation=None,
                 self_loop=False,
                 dropout=0.0):
        super(RelGraphConv, self).__init__()
        self.in_feat = in_feat
        self.out_feat = out_feat
        self.num_rels = num_rels
        self.regularizer = regularizer
        self.num_bases = num_bases
        if self.num_bases is None or self.num_bases > self.num_rels or self.num_bases < 0:
            self.num_bases = self.num_rels
        self.bias = bias
        self.activation = activation
        self.self_loop = self_loop

        if regularizer == "basis":
            # add basis weights
            self.weight = self.params.get(
                'weight', shape=(self.num_bases, self.in_feat, self.out_feat),
                init=mx.init.Xavier(magnitude=math.sqrt(2.0)))
            if self.num_bases < self.num_rels:
                # linear combination coefficients
                self.w_comp = self.params.get(
                    'w_comp', shape=(self.num_rels, self.num_bases),
                    init=mx.init.Xavier(magnitude=math.sqrt(2.0)))
            # message func
            self.message_func = self.basis_message_func
        elif regularizer == "bdd":
            if in_feat % num_bases != 0 or out_feat % num_bases != 0:
                raise ValueError('Feature size must be a multiplier of num_bases.')
            # add block diagonal weights
            self.submat_in = in_feat // self.num_bases
            self.submat_out = out_feat // self.num_bases

            # assuming in_feat and out_feat are both divisible by num_bases
            self.weight = self.params.get(
                'weight',
                shape=(self.num_rels, self.num_bases * self.submat_in * self.submat_out),
                init=mx.init.Xavier(magnitude=math.sqrt(2.0)))
            # message func
            self.message_func = self.bdd_message_func
        else:
            raise ValueError("Regularizer must be either 'basis' or 'bdd'")

        # bias
        if self.bias:
            self.h_bias = self.params.get('bias', shape=(out_feat,),
                                          init=mx.init.Zero())

        # weight for self loop
        if self.self_loop:
            self.loop_weight = self.params.get(
                'W_0', shape=(in_feat, out_feat),
                init=mx.init.Xavier(magnitude=math.sqrt(2.0)))

        self.dropout = nn.Dropout(dropout)

    def basis_message_func(self, edges):
        """Message function for basis regularizer"""
        ctx = edges.src['h'].context
        if self.num_bases < self.num_rels:
            # generate all weights from bases
            weight = self.weight.data(ctx).reshape(
                self.num_bases, self.in_feat * self.out_feat)
            weight = nd.dot(self.w_comp.data(ctx), weight).reshape(
                self.num_rels, self.in_feat, self.out_feat)
        else:
            weight = self.weight.data(ctx)

        msg = utils.bmm_maybe_select(edges.src['h'], weight, edges.data['type'])
        if 'norm' in edges.data:
            msg = msg * edges.data['norm']
        return {'msg': msg}

    def bdd_message_func(self, edges):
        """Message function for block-diagonal-decomposition regularizer"""
        ctx = edges.src['h'].context
        if edges.src['h'].dtype in (np.int32, np.int64) and len(edges.src['h'].shape) == 1:
            raise TypeError('Block decomposition does not allow integer ID feature.')
        weight = self.weight.data(ctx)[edges.data['type'], :].reshape(
            -1, self.submat_in, self.submat_out)
        node = edges.src['h'].reshape(-1, 1, self.submat_in)
        msg = nd.batch_dot(node, weight).reshape(-1, self.out_feat)
        if 'norm' in edges.data:
            msg = msg * edges.data['norm']
        return {'msg': msg}

    def forward(self, g, x, etypes, norm=None):
        r"""Forward computation

        Parameters
        ----------
        g : DGLGraph
            The graph.
        x : mx.ndarray.NDArray
            Input node features. Could be either
              - :math:`(|V|, D)` dense tensor
              - :math:`(|V|,)` int64 vector, representing the categorical values of each
                node. We then treat the input feature as an one-hot encoding feature.
        etypes : mx.ndarray.NDArray
            Edge type tensor. Shape: :math:`(|E|,)`
        norm : mx.ndarray.NDArray
            Optional edge normalizer tensor. Shape: :math:`(|E|, 1)`

        Returns
        -------
        mx.ndarray.NDArray
            New node features.
        """
        g = g.local_var()
        g.ndata['h'] = x
        g.edata['type'] = etypes
        if norm is not None:
            g.edata['norm'] = norm
        if self.self_loop:
            loop_message = utils.matmul_maybe_select(x, self.loop_weight.data(x.context))

        # message passing
        g.update_all(self.message_func, fn.sum(msg='msg', out='h'))

        # apply bias and activation
        node_repr = g.ndata['h']
        if self.bias:
            node_repr = node_repr + self.h_bias.data(x.context)
        if self.self_loop:
            node_repr = node_repr + loop_message
        if self.activation:
            node_repr = self.activation(node_repr)
        node_repr = self.dropout(node_repr)

        return node_repr


text:
# Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*- coding: utf-8 -*-

import pytest
import numpy as np

asdf = pytest.importorskip('asdf', minversion='2.0.0.dev0')
from asdf import util
from asdf.tests import helpers
import astropy.units as u
import astropy.coordinates as coords
from ..fixed_location import FixedLocation
from ..telescope import Telescope
from .test_fixed_location import test1 as fl_test1

def test1(tmpdir, ret=False):

    loc = fl_test1(None, ret=True)
    tel = Telescope(name='VLA', location=loc)
    tree = {'telescope': tel}
    if ret:
        return tel
    helpers.assert_roundtrip_tree(tree, tmpdir)

def test2(tmpdir, ret=False):

    loc = fl_test1(None, ret=True)
    tel = Telescope(name='VLA', location=loc,
                    location_name="Plains of Saint Augustine",
                    telescope_type=['Radio'],
                    org_url='nrao.edu',
                    telescope_url='https://public.nrao.edu/telescopes/vla/',
                    meta={'fines':'hitting cows'})
    tree = {'telescope': tel}
    if ret:
        return tel
    helpers.assert_roundtrip_tree(tree, tmpdir)

ext: py

reason: import numpy,import astropy


text:
from typing import Optional
from zquantum.core.typing import Specs, Union
import numpy as np
from zquantum.core.utils import load_from_specs
from zquantum.core.serialization import save_array
from zquantum.optimizers.grid_search import (
    build_uniform_param_grid as _build_uniform_param_grid,
    save_parameter_grid,
    load_parameter_grid,
    ParameterGrid,
)

# Build uniform parameter grid
def build_uniform_param_grid(
    ansatz_specs: Optional[Specs] = None,
    number_of_params_per_layer: Optional[int] = None,
    number_of_layers: int = 1,
    min_value: float = 0,
    max_value: float = 2 * np.pi,
    step: float = np.pi / 5,
):
    assert (ansatz_specs is None) != (number_of_params_per_layer is None)

    if ansatz_specs is not None:
        ansatz = load_from_specs(ansatz_specs)
        number_of_params = ansatz.number_of_params
    else:
        number_of_params = number_of_params_per_layer

    grid = _build_uniform_param_grid(
        number_of_params, number_of_layers, min_value, max_value, step
    )
    save_parameter_grid(grid, "parameter-grid.json")

def get_parameter_values_list_from_grid(grid: Union[str, ParameterGrid]):
    if isinstance(grid, str):
        grid = load_parameter_grid(grid)

    parameter_values_list = np.array(grid.params_list)

    save_array(parameter_values_list, "parameter-values-list.json")

max_stars_repo_name: zapatacomputing/z-quantum-optimizers

max_issues_repo_name: zapatacomputing/z-quantum-optimizers

max_forks_repo_name: zapatacomputing/z-quantum-optimizers

reason: import numpy


text:
"""
__author__      = "Shivvrat Arya"
__version__     = "Python3.7"
"""
import sys
from copy import deepcopy
from itertools import product as iter_product

from more_itertools import locate
from numpy import prod, product, random, put, sum, divide, array, log10

# todo check if both the type complete and incomplete examples_with_weights give same format
def get_index_given_truth_values(variables, truth_values, cardinalities):
	"""
	The function converts truth table tuples to array indices
	:param variables: The variable in factor
	:param truth_values: The values of given variables in the truth table (same order as variable)
	:param cardinalities: The cardinality of the given variables (same order as variable)
	:return:The index for the array
	"""
	index = 0

	number = 0
	while number < len(variables):
		number_of_tuples_for_current_var = prod(cardinalities[number + 1:]) * truth_values[number]
		index = index + number_of_tuples_for_current_var
		number += 1
	return int(index)

def get_truth_values_given_index(variables, index_value, cardinalities):
	"""
	Gives the truth value (values in tuple) for given index of the array
	@param variables: list containing all the variables in the graph which are in the factor
	@param index_value: Index value for which the tuple is to be founded
	@param cardinalities: Cardinalities of the given variables
	@return: the tuple for corresponding index value
	"""
	number = 0
	truth_table_value = []
	while number < len(variables):
		truth_table_value.append(int(index_value // prod(cardinalities[number + 1:])))
		index_value = index_value - (index_value // prod(cardinalities[number + 1:])) * prod(
				cardinalities[number + 1:])
		number += 1
	return truth_table_value

def generate_random_parameters(var_in_clique, cardinalities_of_var):
	parameters = []
	cardinalities_of_var = array(cardinalities_of_var)
	var_in_clique_np = array(var_in_clique)
	for each_clique in range(len(var_in_clique)):
		var_in_this_clique = var_in_clique_np[each_clique]
		cardinalities_of_var_in_this_clique = cardinalities_of_var[var_in_this_clique]
		total_tuples = product(cardinalities_of_var_in_this_clique)
		if len(var_in_this_clique) > 1:
			# indexing of parameters is done on the basis of parent variables only
			number_of_parameters = int(total_tuples / cardinalities_of_var_in_this_clique[-1])
		else:
			number_of_parameters = 1
		parameters_for_this_clique = random.rand(number_of_parameters)
		if len(parameters_for_this_clique) == 1:
			parameters.append(float(parameters_for_this_clique))
		else:
			parameters.append(parameters_for_this_clique.tolist())
	return parameters

def complete_data(examples, cardinalities_of_var, parameters, var_in_clique):
	final_examples = []
	final_weights = []
	for each_example in examples:
		if -1 not in each_example:
			weight = 1
			final_examples.append(each_example)
			final_weights.append(weight)
		elif -1 in each_example:
			weights_for_completed_examples = []
			completed_examples = []

			var_not_present = list(locate(each_example, lambda a: a == -1))
			cardinalities_of_var = array(cardinalities_of_var)
			cardinalities_of_var_not_present = cardinalities_of_var[var_not_present]
			all_possible_completions = []
			for each_var in range(len(var_not_present)):
				completions_for_var_not_present = list(range(cardinalities_of_var_not_present[each_var]))
				all_possible_completions.append(completions_for_var_not_present)
			all_completions_only_missing_var = list(iter_product(*all_possible_completions))
			for each_tuple_missing_var in all_completions_only_missing_var:
				temp_example = deepcopy(each_example)
				put(temp_example, var_not_present, each_tuple_missing_var)
				weight_for_current_example = 1
				for each_clique in range(len(var_in_clique)):
					if compare_list(var_not_present, var_in_clique[each_clique]):
						values_for_current_clique = temp_example[var_in_clique[each_clique]]
						if len(var_in_clique[each_clique]) > 1:
							index_for_current_example = get_index_given_truth_values(var_in_clique[each_clique][:-1],
							                                                         values_for_current_clique[:-1],
							                                                         cardinalities_of_var[
								                                                         var_in_clique[each_clique][
								                                                         :-1]])
							if values_for_current_clique[-1] == 0:
								weight_for_current_example *= parameters[each_clique][index_for_current_example]
							elif values_for_current_clique[-1] == 1:
								weight_for_current_example *= (1 - parameters[each_clique][index_for_current_example])
							else:
								print("Error")
						elif len(var_in_clique[each_clique]) == 1:
							if values_for_current_clique[0] == 0:
								weight_for_current_example *= parameters[each_clique]
							elif values_for_current_clique[0] == 1:
								weight_for_current_example *= (1 - parameters[each_clique])
							else:
								print("Error")
				weights_for_completed_examples.append(weight_for_current_example)
				completed_examples.append(temp_example)
			denom = sum(weights_for_completed_examples)
			weights_for_completed_examples = divide(weights_for_completed_examples, denom)
			for weight, example in zip(weights_for_completed_examples, completed_examples):
				final_examples.append(example)
				final_weights.append(float(weight))
	return final_examples, final_weights

def compare_list(list_1, list_2):
	for each in list_2:
		if each in list_1:
			return True
	return False

def log_sum_exp(a):
	"""
	Uses the log sum exp technique to do computation in higher domain than float
	@param a: The array for which we want to find the sum
	@return: The sum of the array founded by log sum exp technique
	"""
	# return logsum1(a)/log(10)
	mx = max(a)
	to_return = log10(sum([threshold(10 ** (i - mx)) for i in a])) + mx
	return to_return

def threshold(x):
	"""
	This function is used to find the value in the given limit
	@param x: The value which needs to be checked
	@return: the value according to the threshold
	"""
	if x == float('inf'):
		return sys.float_info.max
	return x

max_stars_repo_name: Shivvrat/Machine-Learning-Algorithms

max_issues_repo_name: Shivvrat/Machine-Learning-Algorithms

max_forks_repo_name: Shivvrat/Machine-Learning-Algorithms

reason: from numpy


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

# Copyright (c) 2016-2021 by University of Kassel and Fraunhofer Institute for Energy Economics
# and Energy System Technology (IEE), Kassel. All rights reserved.

import numpy as np
import pandas as pd

from pandapower.auxiliary import _sum_by_group
from pandapower.pypower.idx_bus import BASE_KV
from pandapower.pypower.idx_gen import GEN_BUS, MBASE
from pandapower.shortcircuit.idx_brch import IKSS_F, IKSS_T, IP_F, IP_T, ITH_F, ITH_T
from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX, KAPPA, R_EQUIV, IKSS1, IP, ITH, X_EQUIV, IKSS2, IKCV, M
from pandapower.shortcircuit.impedance import _calc_zbus_diag

from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci

def _calc_ikss(net, ppc, bus=None):
    # Vectorized for multiple bus
    if bus is None:
        # Slice(None) is equal to : select
        bus_idx = slice(None)
    else:
        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)

    fault = net._options["fault"]
    case = net._options["case"]
    c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
    ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]

    z_equiv = abs(ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j)
    if fault == "3ph":
        ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / np.sqrt(3) * ppc["baseMVA"]
    elif fault == "2ph":
        ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / 2 * ppc["baseMVA"]
    _current_source_current(net, ppc)

def _calc_ikss_1ph(net, ppc, ppc_0, bus=None):
    # Vectorized for multiple bus
    if bus is None:
        # Slice(None) is equal to : select
        bus_idx = slice(None)
    else:
        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)

    case = net._options["case"]
    c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
    ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
    ppc_0["internal"]["baseI"] = ppc_0["bus"][:, BASE_KV] * np.sqrt(3) / ppc_0["baseMVA"]

    z_equiv = abs((ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j) * 2 +
                  (ppc_0["bus"][bus_idx, R_EQUIV] + ppc_0["bus"][bus_idx, X_EQUIV] * 1j))

    ppc_0["bus"][bus_idx, IKSS1] = c / z_equiv / ppc_0["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
    ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppc["baseMVA"]
    _current_source_current(net, ppc)

def _current_source_current(net, ppc):
    ppc["bus"][:, IKCV] = 0
    ppc["bus"][:, IKSS2] = 0
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
    if not False in net.sgen.current_source.values:
        sgen = net.sgen[net._is_elements["sgen"]]
    else:
        sgen = net.sgen[net._is_elements["sgen"] & net.sgen.current_source]
    if len(sgen) == 0:
        return
    if any(pd.isnull(sgen.sn_mva)):
        raise ValueError("sn_mva needs to be specified for all sgens in net.sgen.sn_mva")
    baseI = ppc["internal"]["baseI"]
    sgen_buses = sgen.bus.values
    sgen_buses_ppc = bus_lookup[sgen_buses]
    if not "k" in sgen:
        raise ValueError("Nominal to short-circuit current has to specified in net.sgen.k")
    i_sgen_pu = sgen.sn_mva.values / net.sn_mva * sgen.k.values
    buses, ikcv_pu, _ = _sum_by_group(sgen_buses_ppc, i_sgen_pu, i_sgen_pu)
    ppc["bus"][buses, IKCV] = ikcv_pu
    if net["_options"]["inverse_y"]:
        Zbus = ppc["internal"]["Zbus"]
        ppc["bus"][:, IKSS2] = abs(1 / np.diag(Zbus) * np.dot(Zbus, ppc["bus"][:, IKCV] * -1j) / baseI)
    else:
        ybus_fact = ppc["internal"]["ybus_fact"]
        diagZ = _calc_zbus_diag(net, ppc)
        ppc["bus"][:, IKSS2] = abs(ybus_fact(ppc["bus"][:, IKCV] * -1j) / diagZ / baseI)
    ppc["bus"][buses, IKCV] /= baseI[buses]

def _calc_ip(net, ppc):
    ip = np.sqrt(2) * (ppc["bus"][:, KAPPA] * ppc["bus"][:, IKSS1] + ppc["bus"][:, IKSS2])
    ppc["bus"][:, IP] = ip

def _calc_ith(net, ppc):
    tk_s = net["_options"]["tk_s"]
    kappa = ppc["bus"][:, KAPPA]
    f = 50
    n = 1
    m = (np.exp(4 * f * tk_s * np.log(kappa - 1)) - 1) / (2 * f * tk_s * np.log(kappa - 1))
    m[np.where(kappa > 1.99)] = 0
    ppc["bus"][:, M] = m
    ith = (ppc["bus"][:, IKSS1] + ppc["bus"][:, IKSS2]) * np.sqrt(m + n)
    ppc["bus"][:, ITH] = ith

def _calc_ib_generator(net, ppci):
    Zbus = ppci["internal"]["Zbus"]
    baseI = ppci["internal"]["baseI"]
    tk_s = net._options['tk_s']
    c = 1.1

    z_equiv = ppci["bus"][:, R_EQUIV] + ppci["bus"][:, X_EQUIV] * 1j
    I_ikss = c / z_equiv / ppci["bus"][:, BASE_KV] / np.sqrt(3) * ppci["baseMVA"]

    # calculate voltage source branch current
    # I_ikss = ppci["bus"][:, IKSS1]
    V_ikss = (I_ikss * baseI) * Zbus

    gen = net["gen"][net._is_elements["gen"]]
    gen_vn_kv = gen.vn_kv.values

    gen_buses = ppci['gen'][:, GEN_BUS].astype(np.int64)
    gen_mbase = ppci['gen'][:, MBASE]
    gen_i_rg = gen_mbase / (np.sqrt(3) * gen_vn_kv)

    gen_buses_ppc, gen_sn_mva, I_rG = _sum_by_group(gen_buses, gen_mbase, gen_i_rg)

    # shunt admittance of generator buses and generator short circuit current
    # YS = ppci["bus"][gen_buses_ppc, GS] + ppci["bus"][gen_buses_ppc, BS] * 1j
    # I_kG = V_ikss.T[:, gen_buses_ppc] * YS / baseI[gen_buses_ppc]

    xdss_pu = gen.xdss_pu.values
    rdss_pu = gen.rdss_pu.values
    cosphi = gen.cos_phi.values
    X_dsss = xdss_pu * np.square(gen_vn_kv) / gen_mbase
    R_dsss = rdss_pu * np.square(gen_vn_kv) / gen_mbase

    K_G = ppci['bus'][gen_buses, BASE_KV] / gen_vn_kv * c / (1 + xdss_pu * np.sin(np.arccos(cosphi)))
    Z_G = (R_dsss + 1j * X_dsss)

    I_kG = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3) / (Z_G * K_G) * ppci["baseMVA"]

    dV_G = 1j * X_dsss * K_G * I_kG
    V_Is = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3)

    # I_kG_contribution = I_kG.sum(axis=1)
    # ratio_SG_ikss = I_kG_contribution / I_ikss
    # close_to_SG = ratio_SG_ikss > 5e-2

    close_to_SG = I_kG / I_rG > 2

    if tk_s == 2e-2:
        mu = 0.84 + 0.26 * np.exp(-0.26 * abs(I_kG) / I_rG)
    elif tk_s == 5e-2:
        mu = 0.71 + 0.51 * np.exp(-0.3 * abs(I_kG) / I_rG)
    elif tk_s == 10e-2:
        mu = 0.62 + 0.72 * np.exp(-0.32 * abs(I_kG) / I_rG)
    elif tk_s >= 25e-2:
        mu = 0.56 + 0.94 * np.exp(-0.38 * abs(I_kG) / I_rG)
    else:
        raise UserWarning('not implemented for other tk_s than 20ms, 50ms, 100ms and >=250ms')

    mu = np.clip(mu, 0, 1)

    I_ikss_G = abs(I_ikss - np.sum((1 - mu) * I_kG, axis=1))

    # I_ikss_G = I_ikss - np.sum(abs(V_ikss.T[:, gen_buses_ppc]) * (1-mu) * I_kG, axis=1)

    I_ikss_G = abs(I_ikss - np.sum(dV_G / V_Is * (1 - mu) * I_kG, axis=1))

    return I_ikss_G

def _calc_single_bus_sc(net, ppc, bus):
    # case = net._options["case"]
    bus_idx = net._pd2ppc_lookups["bus"][bus]
    n = ppc["bus"].shape[0]
    Zbus = ppc["internal"]["Zbus"]
    #    Yf = ppc["internal"]["Yf"]
    #    Yt = ppc["internal"]["Yf"]
    baseI = ppc["internal"]["baseI"]
    #    fb = np.real(ppc["branch"][:, 0]).astype(int)
    #    tb = np.real(ppc["branch"][:, 1]).astype(int)
    # c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]

    # calculate voltage source branch current
    V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
    V = V_ikss[:, bus_idx]
    #    ikss_all_f = np.conj(Yf.dot(V_ikss))
    #    ikss_all_t = np.conj(Yt.dot(V_ikss))
    current_sources = any(ppc["bus"][:, IKCV]) > 0
    if current_sources:
        current = np.tile(-ppc["bus"][:, IKCV], (n, 1))
        np.fill_diagonal(current, current.diagonal() + ppc["bus"][:, IKSS2])
        V_source = np.dot((current * baseI), Zbus).T
        V = V + V_source[:, bus_idx]
    # add current source branch current if there is one
    #    ppc["branch"][:, IKSS_F] = abs(ikss_all_f[:, bus_idx] / baseI[fb])
    #    ppc["branch"][:, IKSS_T] = abs(ikss_all_t[:, bus_idx] / baseI[tb])
    calc_branch_results(net, ppc, V)

def _calc_single_bus_sc_no_y_inv(net, ppc, bus):
    # Vectorized for multiple bus
    if bus is None:
        # Slice(None) is equal to : select
        bus_idx = slice(None)
    else:
        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)

    ybus = ppc["internal"]["Ybus"]
    ybus_fact = ppc["internal"]["ybus_fact"]
    # case = net._options["case"]
    baseI = ppc["internal"]["baseI"]
    # vqj = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX] #this is the source voltage in per unit (VQj)

    # Solve Ikss from voltage source
    n_bus = ybus.shape[0]

    # ybus_sub_mask = (np.arange(ybus.shape[0]) != bus_idx)
    # V_ikss = np.zeros(n_bus, dtype=np.complex)
    # V_ikss[bus_idx] = vqj[bus_idx]

    # # Solve Ax = b
    # b = np.zeros(n_bus-1, dtype=np.complex) -\
    #     (ybus[:, ~ybus_sub_mask].toarray())[ybus_sub_mask].ravel() * V_ikss[bus_idx]
    # ybus_sub = ybus[ybus_sub_mask, :][:, ybus_sub_mask]
    # x = spsolve(ybus_sub, b)

    # V_ikss[ybus_sub_mask] = x
    # I_ikss = np.zeros(n_bus, dtype=np.complex)
    # I_ikss[bus_idx] = np.dot(ybus[bus_idx, :].toarray(), V_ikss)
    # V = V_ikss

    # Version 2
    I_ikss = np.zeros(n_bus, dtype=np.complex)
    I_ikss[bus_idx] = ppc["bus"][bus_idx, IKSS1]
    V_ikss = ybus_fact(I_ikss * baseI)
    V = V_ikss

    #TODO include current sources
    current_sources = any(ppc["bus"][:, IKCV]) > 0
    if current_sources:
        current = -ppc["bus"][:, IKCV]
        current[bus_idx] += ppc["bus"][bus_idx, IKSS2]
        V_source = ybus_fact(current)
        V += V_source

    calc_branch_results(net, ppc, V)

def calc_branch_results(net, ppci, V):
    Ybus = ppci["internal"]["Ybus"]
    Yf = ppci["internal"]["Yf"]
    Yt = ppci["internal"]["Yt"]
    baseMVA, bus, gen, branch, ref, _, _, _, _, _, ref_gens = _get_pf_variables_from_ppci(ppci)
    bus, gen, branch = pfsoln_pypower(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens)
    ppci["bus"], ppci["gen"], ppci["branch"] = bus, gen, branch

def _calc_branch_currents(net, ppc, bus):
    # Vectorized for multiple bus
    if bus is None:
        # Slice(None) is equal to select all
        bus = net.bus.index

    bus_idx = net._pd2ppc_lookups["bus"][bus]
    # Select only in service bus for sc calculation
    bus_idx = bus_idx[bus_idx < ppc['bus'].shape[0]]
    n_sc_bus = np.shape(bus_idx)[0]

    case = net._options["case"]

    Yf = ppc["internal"]["Yf"]
    Yt = ppc["internal"]["Yt"]
    baseI = ppc["internal"]["baseI"]
    n_bus = ppc["bus"].shape[0]
    fb = np.real(ppc["branch"][:, 0]).astype(int)
    tb = np.real(ppc["branch"][:, 1]).astype(int)
    minmax = np.nanmin if case == "min" else np.nanmax

    # calculate voltage source branch current
    if net["_options"]["inverse_y"]:
        Zbus = ppc["internal"]["Zbus"]
        V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
        V_ikss = V_ikss[:, bus_idx]
    else:
        ybus_fact = ppc["internal"]["ybus_fact"]
        V_ikss = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
        for ix, b in enumerate(bus_idx):
            ikss = np.zeros(n_bus, dtype=np.complex)
            ikss[b] = ppc["bus"][b, IKSS1] * baseI[b]
            V_ikss[:, ix] = ybus_fact(ikss)

    ikss1_all_f = np.conj(Yf.dot(V_ikss))
    ikss1_all_t = np.conj(Yt.dot(V_ikss))
    ikss1_all_f[abs(ikss1_all_f) < 1e-10] = 0.
    ikss1_all_t[abs(ikss1_all_t) < 1e-10] = 0.

    # add current source branch current if there is one
    current_sources = any(ppc["bus"][:, IKCV]) > 0
    if current_sources:
        current = np.tile(-ppc["bus"][:, IKCV], (n_sc_bus, 1))
        for ix, b in enumerate(bus_idx):
            current[ix, b] += ppc["bus"][b, IKSS2]

        # calculate voltage source branch current
        if net["_options"]["inverse_y"]:
            Zbus = ppc["internal"]["Zbus"]
            V = np.dot((current * baseI), Zbus).T
        else:
            ybus_fact = ppc["internal"]["ybus_fact"]
            V = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
            for ix, b in enumerate(bus_idx):
                V[:, ix] = ybus_fact(current[ix, :] * baseI[b])

        fb = np.real(ppc["branch"][:, 0]).astype(int)
        tb = np.real(ppc["branch"][:, 1]).astype(int)
        ikss2_all_f = np.conj(Yf.dot(V))
        ikss2_all_t = np.conj(Yt.dot(V))

        ikss_all_f = abs(ikss1_all_f + ikss2_all_f)
        ikss_all_t = abs(ikss1_all_t + ikss2_all_t)
    else:
        ikss_all_f = abs(ikss1_all_f)
        ikss_all_t = abs(ikss1_all_t)

    if net._options["return_all_currents"]:
        ppc["internal"]["branch_ikss_f"] = ikss_all_f / baseI[fb, None]
        ppc["internal"]["branch_ikss_t"] = ikss_all_t / baseI[tb, None]
    else:
        ikss_all_f[abs(ikss_all_f) < 1e-10] = np.nan
        ikss_all_t[abs(ikss_all_t) < 1e-10] = np.nan
        ppc["branch"][:, IKSS_F] = np.nan_to_num(minmax(ikss_all_f, axis=1) / baseI[fb])
        ppc["branch"][:, IKSS_T] = np.nan_to_num(minmax(ikss_all_t, axis=1) / baseI[tb])

    if net._options["ip"]:
        kappa = ppc["bus"][:, KAPPA]
        if current_sources:
            ip_all_f = np.sqrt(2) * (ikss1_all_f * kappa[bus_idx] + ikss2_all_f)
            ip_all_t = np.sqrt(2) * (ikss1_all_t * kappa[bus_idx] + ikss2_all_t)
        else:
            ip_all_f = np.sqrt(2) * ikss1_all_f * kappa[bus_idx]
            ip_all_t = np.sqrt(2) * ikss1_all_t * kappa[bus_idx]

        if net._options["return_all_currents"]:
            ppc["internal"]["branch_ip_f"] = abs(ip_all_f) / baseI[fb, None]
            ppc["internal"]["branch_ip_t"] = abs(ip_all_t) / baseI[tb, None]
        else:
            ip_all_f[abs(ip_all_f) < 1e-10] = np.nan
            ip_all_t[abs(ip_all_t) < 1e-10] = np.nan
            ppc["branch"][:, IP_F] = np.nan_to_num(minmax(abs(ip_all_f), axis=1) / baseI[fb])
            ppc["branch"][:, IP_T] = np.nan_to_num(minmax(abs(ip_all_t), axis=1) / baseI[tb])

    if net._options["ith"]:
        n = 1
        m = ppc["bus"][bus_idx, M]
        ith_all_f = ikss_all_f * np.sqrt(m + n)
        ith_all_t = ikss_all_t * np.sqrt(m + n)

        if net._options["return_all_currents"]:
            ppc["internal"]["branch_ith_f"] = ith_all_f / baseI[fb, None]
            ppc["internal"]["branch_ith_t"] = ith_all_t / baseI[tb, None]
        else:
            ppc["branch"][:, ITH_F] = np.nan_to_num(minmax(ith_all_f, axis=1) / baseI[fb])
            ppc["branch"][:, ITH_T] = np.nan_to_num(minmax(ith_all_t, axis=1) / baseI[fb])

max_stars_repo_name: julffers/pandapower

max_issues_repo_name: julffers/pandapower

max_forks_repo_name: julffers/pandapower


text:
#coding: utf-8
import argparse
import numpy as np
import chainer
from chainer import cuda
import chainer.functions as F
from chainer import optimizers
import time
import pickle

import random
import os
import cv2

parser = argparse.ArgumentParser(description='Chainer example: MNIST')
parser.add_argument('--gpu', '-gpu', default=-1, type=int,
                    help='GPU ID (negative value indicates CPU)')

# GPUが使えるか確認
args = parser.parse_args()
if args.gpu >= 0:
    cuda.check_cuda_available()
xp = cuda.cupy if args.gpu >= 0 else np

# 学習のパラメータ
batchsize = 1
n_epoch = 1
n_units = 512
n_input = 100 * 100 * 3

# テスト精度を出力するファイル名
accfile = "accuracy.txt"

# データセットをロード
print 'load BABY dataset'
t_data = np.load("t_data.npy").astype(np.float32)
f_data = np.load("f_data.npy").astype(np.float32)
t_label = np.ones((t_data.shape[0])).astype(np.int32)
f_label = np.zeros((f_data.shape[0])).astype(np.int32)

N_all_data = t_data.shape[0] + f_data.shape[0]

x_train_a, x_test_a = np.split(t_data, [int(t_data.shape[0] * 0.8)])
x_train_b, x_test_b = np.split(f_data, [int(f_data.shape[0] * 0.8)])

x_train = np.r_[x_train_a, x_train_b]
x_test = np.r_[x_test_a, x_test_b]

y_train_a, y_test_a = np.split(t_label, [int(t_data.shape[0] * 0.8)])
y_train_b, y_test_b = np.split(f_label, [int(f_data.shape[0] * 0.8)])
y_train = np.r_[y_train_a, y_train_b]
y_test = np.r_[y_test_a, y_test_b]

N = y_train.size
N_test = y_test.size

print N, N_test, N_all_data

# 多層パーセプトロンのモデル（パラメータ集合）
model = pickle.load(open("model", "rb"))

# GPU使用のときはGPUにモデルを転送
if args.gpu >= 0:
    cuda.get_device(args.gpu).use()
    model.to_gpu()

def forward(x_data, y_data, train=True):
    # 順伝播の処理を定義
    x = chainer.Variable(x_data)
    # 隠れ層1の出力
    h1 = F.dropout(F.relu(model.l1(x)), train=train)
    # 隠れ層2の出力
    h2 = F.dropout(F.relu(model.l2(h1)), train=train)
    # 出力層の出力
    y = model.l3(h2)

    return y

# Optimizerをセット
# 最適化対象であるパラメータ集合のmodelを渡しておく
optimizer = optimizers.Adam()
optimizer.setup(model)

src_dir = "image/"
files = os.listdir(src_dir)
files.remove(".DS_Store")

i = 0
for i in xrange(1000):

    file = random.choice(files)
    print file
    img = cv2.imread(src_dir+file)
    resized_x = 512
    resized_y = int(512 / img.shape[0] * img.shape[1])

    resized_img = cv2.resize(img, (resized_x, resized_y))
    patch_size = 128
    x = random.randint(0, resized_x-patch_size)
    y = random.randint(0, resized_y-patch_size)
    data = resized_img[y:y+patch_size, x:x+patch_size]

    label = np.argmax(forward(np.array([data.flatten()], dtype=np.float32),
                              None, train=False).data)

    if label == 1:
        #img = np.array(x_train[i].reshape(100, 100, 3))*255
        #img = img.astype(np.int32)
        cv2.imwrite("detected_unchi/t/"+str(i)+".bmp", data)
    elif label == 0:
        #img = np.array(x_train[i].reshape(100, 100, 3))*255
        #img = img.astype(np.int32)
        cv2.imwrite("detected_unchi/f/"+str(i)+".bmp", data)

max_stars_repo_name: wbawakate/baby-unchi

max_issues_repo_name: wbawakate/baby-unchi

max_forks_repo_name: wbawakate/baby-unchi


text:
import copy
import random

from player import Player
import numpy as np
import json

class Evolution:
    def __init__(self):
        self.game_mode = "Neuroevolution"
        fitness_results = {
                'min_fitness': [],
                'max_fitness': [],
                'mean_fitness': []
            }
        with open('fitness_results.json', 'w') as out_file:
                json.dump(fitness_results, out_file)
    def clone_player(self, player):
        """
        Gets a player as an input and produces a clone of that player.
        """
        new_player = Player(self.game_mode)
        new_player.nn = player.nn
        new_player.fitness = player.fitness
        return new_player

    #By this function we calculate
    def calculate_cumulative_probabilities(self, players):
        total_fitness = 0
        for player in players:
            total_fitness += player.fitness
        probabilities = []
        for player in players:
            probabilities.append(player.fitness / total_fitness)
        # convert it to
        for i in range(1, len(players)):
            probabilities[i] += probabilities[i - 1]
        return probabilities
    def roulette_wheel(self, players, parent_numbers):
        probabilities = self.calculate_cumulative_probabilities(players)

        results = []
        for random_number in np.random.uniform(low=0, high=1, size=parent_numbers):
            for i, probability in enumerate(probabilities):
                if random_number <= probability:
                    results.append(self.clone_player(players[i]))
                    break

        return results

    def save_fitness_result(self, min_fitness, max_fitness, mean_fitness):
        fitness_results = {
                'min_fitness': [],
                'max_fitness': [],
                'mean_fitness': []
            }
        with open('fitness_results.json', 'r') as in_file:
                fitness_results = json.load(in_file)

        fitness_results['min_fitness'].append(min_fitness)
        fitness_results['max_fitness'].append(max_fitness)
        fitness_results['mean_fitness'].append(mean_fitness)

        with open('fitness_results.json', 'w') as out_file:
                json.dump(fitness_results, out_file)
    def next_population_selection(self, players, num_players):
        """
        Gets list of previous and current players (μ + λ) and returns num_players number of players based on their
        fitness value.

        :param players: list of players in the previous generation
        :param num_players: number of players that we return
        """
        result = players
        selection_method = "SUS"

        if selection_method == "top-k":
            sorted_players = sorted(players, key=lambda player: player.fitness, reverse=True)
            result = sorted_players[: num_players]
        elif selection_method == 'roulette wheel':
            result = self.roulette_wheel(players, num_players)
        elif selection_method == "SUS":
            result = self.sus(players, num_players)

        fitness_list = [player.fitness for player in players]
        max_fitness = float(np.max(fitness_list))
        mean_fitness = float(np.mean(fitness_list))
        min_fitness = float(np.min(fitness_list))
        self.save_fitness_result(min_fitness, max_fitness, mean_fitness)
        return result
    def sus(self, players, num_players):
        # Create Intervals
        interval_length = 1 - 1 / num_players
        intervals = np.linspace(0, interval_length, num_players)
        random_number = np.random.uniform(0, 1 / num_players, 1)
        intervals += random_number

        probabilities = self.calculate_cumulative_probabilities(players)

        result = []
        for interval in intervals:
            for i, probability in enumerate(probabilities):
                if interval < probability:
                    result.append(self.clone_player(players[i]))
                    break

        return result
    def mutate_child(self, child):
        number_of_layer = random.randint(0, len(child.nn.weights)-1)
        number_of_weight = random.randint(0, len(child.nn.weights[number_of_layer])-1)
        child.nn.weights[number_of_layer][number_of_weight] = random.random()
        return child

    def q_tournament(self, players, q):
        q_selected = np.random.choice(players, q)
        return max(q_selected, key=lambda player: player.fitness)

    def crossover(self, child1_array, child2_array, parent1_array, parent2_array):
        row_size, column_size = child1_array.shape
        section_1, section_2, section_3 = int(row_size / 3), int(2 * row_size / 3), row_size

        random_number = np.random.uniform(0, 1, 1)
        if random_number > 0.5:
            child1_array[:section_1, :] = parent1_array[:section_1:, :]
            child1_array[section_1:section_2, :] = parent2_array[section_1:section_2, :]
            child1_array[section_2:, :] = parent1_array[section_2:, :]

            child2_array[:section_1, :] = parent2_array[:section_1:, :]
            child2_array[section_1:section_2, :] = parent1_array[section_1:section_2, :]
            child2_array[section_2:, :] = parent2_array[section_2:, :]
        else:
            child1_array[:section_1, :] = parent2_array[:section_1:, :]
            child1_array[section_1:section_2, :] = parent1_array[section_1:section_2, :]
            child1_array[section_2:, :] = parent2_array[section_2:, :]

            child2_array[:section_1, :] = parent1_array[:section_1:, :]
            child2_array[section_1:section_2, :] = parent2_array[section_1:section_2, :]
            child2_array[section_2:, :] = parent1_array[section_2:, :]

    def add_gaussian_noise(self, array, threshold):
        random_number = np.random.uniform(0, 1, 1)
        if random_number < threshold:
            array += np.random.normal(size=array.shape)
    def mutate(self, child):
        # child: an object of class `Player`
        threshold = 0.2
        self.add_gaussian_noise(child.nn.W1, threshold)
        self.add_gaussian_noise(child.nn.W2, threshold)
        self.add_gaussian_noise(child.nn.b1, threshold)
        self.add_gaussian_noise(child.nn.b2, threshold)
    def reproduction(self, parent1, parent2):
        child1 = Player(self.game_mode)
        child2 = Player(self.game_mode)
        self.crossover(child1.nn.W1, child2.nn.W1, parent1.nn.W1, parent2.nn.W1)
        self.crossover(child1.nn.W2, child2.nn.W2, parent1.nn.W2, parent2.nn.W2)
        self.crossover(child1.nn.b1, child2.nn.b1, parent1.nn.b1, parent2.nn.b1)
        self.crossover(child1.nn.b2, child2.nn.b2, parent1.nn.b2, parent2.nn.b2)

        self.mutate(child1)
        self.mutate(child2)
        return child1, child2

    def generate_new_population(self, num_players, prev_players=None):

        if prev_players is None:
            return [Player(self.game_mode) for _ in range(num_players)]

        else:
            method = "Q tournament"
            children = []
            parents = []

            if method == 'roulette wheel':
                parents = self.roulette_wheel(prev_players, num_players)
            elif method == "Q tournament":
                for _ in range(num_players):
                    parents.append(self.q_tournament(prev_players, q=3))

            for i in range(0, len(parents), 2):
                child1, child2 = self.reproduction(parents[i], parents[i + 1])
                children.append(child1)
                children.append(child2)
            return children

max_stars_repo_name: smoeina/Snail-Jumper-Neuroevolution

max_issues_repo_name: smoeina/Snail-Jumper-Neuroevolution

max_forks_repo_name: smoeina/Snail-Jumper-Neuroevolution


text:
"""
    Code by Tae-Hwan Hung(@graykode)
    https://en.wikipedia.org/wiki/Normal_distribution
"""
import numpy as np
from matplotlib import pyplot as plt

def gaussian(x, n):
    u = x.mean()
    s = x.std()

    # divide [x.min(), x.max()] by n
    x = np.linspace(x.min(), x.max(), n)

    a = ((x - u) ** 2) / (2 * (s ** 2))
    y = 1 / (s * np.sqrt(2 * np.pi)) * np.exp(-a)

    return x, y, x.mean(), x.std()

x = np.arange(-100, 100) # define range of x
x, y, u, s = gaussian(x, 10000)

plt.plot(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f$' % (u, s))
plt.legend()
plt.savefig('graph/gaussian.png')
plt.show()


text:
# Copyright 2021 The Distla Authors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# =============================================================================
from jax import lax
import jax.numpy as jnp
import numpy as np

from distla_core.analysis.benchmarks import benchmark_summa
from distla_core.utils import pops

if __name__ == "__main__":
  local_rows = np.array([128, 256, 384, 512, 768, 1024, 1536, 2048, 4096, 8192,
                         16384])
  OUT_ROWS = tuple(pops.NROWS * local_rows)
  OUT_COLS = (None,)  # columns of the output matrix (=out_rows if None)
  SHARED_DIM = (None,)  # shared dim of the multiply (=out_rows if None)
  DTYPES = (jnp.float32,)  # dtype of the matrices
  P_SZS = (128, 256, -1)  # panel size of the multiplication
  TRANSPOSE_AS = (False, True) # whether to do C = A.T @ B
  TRANSPOSE_BS = (False,) # whether to do C = A @ B.T
  PRECISIONS = (lax.Precision.DEFAULT, lax.Precision.HIGH,
                lax.Precision.HIGHEST,) # ASIC matmul precision
  SEED_AS = (None,)  # Random seed to initialize A; system clock if None.
  SEED_BS = (None,)  # Same to initialize B; system clock + 1 if None.
  BATCH_SIZE = 3
  REDUCTION_MODE = "min"  # how to assemblage results
  OUTPUT_DIR_PATH = None  # directory of output; CWD if None
  OUTPUT_NAME = "benchmark_summa" # output saved to OUTPUT_NAME.csv
  _ = benchmark_summa.benchmark_summa(
    OUT_ROWS, out_cols=OUT_COLS, shared_dim=SHARED_DIM, dtypes=DTYPES,
    p_szs=P_SZS, transpose_as=TRANSPOSE_AS, transpose_bs=TRANSPOSE_BS,
    precisions=PRECISIONS, seed_as=SEED_AS, seed_bs=SEED_BS,
    batch_size=BATCH_SIZE, reduction_mode=REDUCTION_MODE,
    output_dir_path=OUTPUT_DIR_PATH, output_name=OUTPUT_NAME)

reason: import numpy,import jax,from jax


text:
#coding:utf-8

##############################################
# iCub learning program
# Basically, Gibbs sampling of GMMs and LDA
# Akira Taniguchi 2016/05/31-2020/07/15
##############################################

#numpy.random.multivariate_normal
#http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.random.multivariate_normal.html
#scipy.stats.multivariate_normal
#http://docs.scipy.org/doc/scipy-0.17.0/reference/generated/scipy.stats.multivariate_normal.html
#scipy.stats.invwishart
#http://docs.scipy.org/doc/scipy-0.17.0/reference/generated/scipy.stats.invwishart.html
#numpy.random.dirichlet
#http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.random.dirichlet.html
#scipy.stats.dirichlet
#http://docs.scipy.org/doc/scipy-0.17.0/reference/generated/scipy.stats.dirichlet.html
#numpy.random.multinomial
#http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.random.multinomial.html
#scipy.stats.rv_discrete
#http://docs.scipy.org/doc/scipy-0.17.0/reference/generated/scipy.stats.rv_discrete.html

import os
import random
import collections
import itertools
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from numpy.random import multinomial,uniform,dirichlet
from scipy.stats import multivariate_normal,invwishart,rv_discrete
from math import pi as PI
from math import cos,sin,sqrt,exp,log,fabs,fsum,degrees,radians,atan2
from sklearn.cluster import KMeans
from __init__ import *

def Makedir(dir):
    try:
        os.mkdir( dir )
    except:
        pass

def stick_breaking(alpha, k):
    betas = np.random.beta(1, alpha, k)
    remaining_pieces = np.append(1, np.cumprod(1 - betas[:-1]))
    p = betas * remaining_pieces
    return p/p.sum()

def Sample_Frame_Random(num):
    #modal = [ modality[i] for i in range(len(modality))]
    F = []
    for i in xrange(num):
      mo = int(uniform(0,len(modality)))
      F = F + [ modality[mo] ]
    return F

def Sample_Frame(num):
    modal = ["a","p","o","c"]
    other = "x"
    dis = 0
    F = []
    if (num >= len(modal)):
      #print "[Error] over modality"
      dis = num - len(modal)
      num = 4
      #return F
    for n in xrange(num):
      mo = int(uniform(0,4-n))
      F = F + [ modal[mo] ]
      modal.pop(mo)
    for i in xrange(dis):
      F.insert( int(uniform( 0, int(num+i+1) )), other )
    return F

def data_read(trialname,finename,sn,en):
    foldername = datafolder + trialname

    M = [0 for d in xrange(D)]   # the set of the number of objects
    N = [0 for d in xrange(D)]   # the set of the number of words
    Ad = [0 for d in xrange(D)]

    w_dn = [ [] for d in xrange(D) ]  # the set of sentences # The sentence is the set of words.
    a_d  = [ [] for d in xrange(D) ]  # the set of actions
    o_dm = [ [] for d in xrange(D) ]  # the set of object features
    c_dm = [ [] for d in xrange(D) ]  # the set of object colors
    p_dm = [ [] for d in xrange(D) ]  # the set of pobject positions

    min_a = [10,10,10]     # Temporary initial value
    max_a = [-10,-10,-10]  # Temporary initial value
    min_o = 10000
    max_o = -10000

    ## reading the set of the number of objects for each sentences
    for d in xrange(D):
      line_num = 0
      ## Read the number of objects M[d]
      for line in open(foldername+str(d+sn).zfill(3)+'/object_center.txt','r'):
        if line_num == 0:
          M[d] = int(line)
          o_dm[d] = [ [ 0 for k in xrange(k_sift) ] for m in xrange(M[d]) ]
          c_dm[d] = [ [ 0 for k in xrange(k_rgb)  ] for m in xrange(M[d]) ]
          p_dm[d] = [ [ 0 for k in xrange(dim_p)  ] for m in xrange(M[d]) ]
          #elif line_num == 1:
          #  Ad[d] = int(line)-1 #物体番号1からMを0からM-1にする
        elif line_num == 1:
          dummy = 0
          #print "if random_obj was changed, but no problem."
        else:
          itemList = line[:-1].split(',')
          #for i in xrange(len(itemList)):
          if itemList[0] != '':
            ## Read the position p_dm
            p_dm[d][int(itemList[0])-1] = [float(itemList[1]),float(itemList[2])]
        line_num = line_num + 1

      gyo = 0
      ## Read the target object position and end-effector position (Hand posiiton of iCub)
      for line in open(foldername+str(d+sn).zfill(3)+'/target_object.txt','r'):
        if gyo == 0:
          Ad[d] = int(line)-1
        elif gyo == 1:
          itemList = line[:-1].split(',')
          obj_pos = np.array([float(itemList[0]),float(itemList[1]),float(itemList[2])])
        elif gyo == 2:
          itemList = line[:-1].split(',')
          enf_pos = np.array([float(itemList[0]),float(itemList[1]),float(itemList[2])])
        elif gyo == 4:
          randomove = float(line)
        gyo = gyo + 1

      tmp = enf_pos - obj_pos #obj_pos - enf_pos
      if (d == 0):
         min_a = list(tmp)
         max_a = list(tmp)
      for i in xrange(3):
        if (min_a[i] > tmp[i]):
           min_a[i]= tmp[i]
        if (max_a[i] < tmp[i]):
           max_a[i] = tmp[i]
      a_d[d] = list(tmp) + [randomove] # Relative 3D position, finger bending (相対3次元位置、指の曲げ具合)

      gyo = 0
      for line in open(foldername+str(d+sn).zfill(3)+'/action.csv','r'):
        itemList = line[:-1].split(',')

        if ("" in itemList):
          itemList.pop(itemList.index(""))
        # Normalized with minimum and maximum values for each joint location (関節箇所ごとに最小値と最大値で正規化)
        if(gyo == 1 or gyo == 7 or gyo == 10):
          for i in xrange(len(itemList)):
            if gyo == 1: #head
              min = [-40,-70,-55,-35,-50, 0]
              max = [ 30, 60, 55, 15, 50,90]
            elif gyo == 7: #right_arm
              min = [-95,    0,-37,15.5,-90,-90,-20,  0,9.6,  0,  0,  0,  0,  0,   0,  0]
              max = [ 10,160.8, 80, 106, 90,  0, 40, 60, 90, 90,180, 90, 180, 90,180,270]
            elif gyo == 10: #torso
              min = [-50,-30,-10]
              max = [ 50, 30, 70]
            if itemList[i] != '':
              #print d,itemList[i]
              a_d[d] = a_d[d] + [ (float(itemList[i])-min[i])/float(max[i]-min[i]) ]  #Normalize and add to array
        if(gyo == 13):
          tactile = [0.0 for i in xrange(len(itemList)/12)]
          for i in xrange(len(itemList)/12):
            #for j in xrange(12):
            #  print i*12+j
            tactile[i] = sum([float(itemList[i*12+j]) for j in xrange(12)])/float(255*12)  #Average of 12 data

        gyo = gyo + 1
      a_d[d] = a_d[d] + tactile

      for m in xrange(M[d]):
        ##物体情報 BoF(SIFT)orCNN特徴の読み込み
        for line in open(foldername+str(d+sn).zfill(3)+'/image/object_' + Descriptor + '_'+str(m+1).zfill(2)+'.csv', 'r'):
          itemList = line[:-1].split(',')
          #print c
          #W_index = W_index + [itemList]
          for i in xrange(len(itemList)):
            if itemList[i] != '':
              #print c,i,itemList[i]
              o_dm[d][m][i] = float(itemList[i])
              if min_o > o_dm[d][m][i]:
                min_o = o_dm[d][m][i]
              if max_o < o_dm[d][m][i]:
                max_o = o_dm[d][m][i]

        if (CNNmode == 0) or (CNNmode == -1):
          # Normalized BoF count (BoFのカウント数を正規化)
          sum_o_dm = sum(o_dm[d][m])
          for i in xrange(len(o_dm[d][m])):
            o_dm[d][m][i] = o_dm[d][m][i] / sum_o_dm

        ## Read color information BoF(RGB)
        for line in open(foldername+str(d+sn).zfill(3)+'/image/object_RGB_BoF_'+str(m+1).zfill(2)+'.csv', 'r'):
          itemList = line[:-1].split(',')
          #print c
          #W_index = W_index + [itemList]
          for i in xrange(len(itemList)):
            if itemList[i] != '':
              #print c,i,itemList[i]
              c_dm[d][m][i] = float(itemList[i])

        # Normalized BoF count (BoFのカウント数を正規化)
        sum_c_dm = sum(c_dm[d][m])
        for i in xrange(len(c_dm[d][m])):
          c_dm[d][m][i] = c_dm[d][m][i] / sum_c_dm

    # Get the word data
    d = 0
    for line in open(foldername +"("+str(sn).zfill(3)+"-"+str(en).zfill(3)+")"+'/' + trialname +'_words.csv', 'r'):
        itemList = line[:-1].split(',')
        for i in xrange(len(itemList)):
            if itemList[i] != '':
              w_dn[d] = w_dn[d] + [str(itemList[i])]
              N[d] = N[d] + 1
        d = d + 1

    print w_dn

    ## Normalize relative coordinates of action data (アクションデータの相対座標を正規化)
    for d in xrange(D):
      a_d[d][0] = (a_d[d][0] - min_a[0]) / (max_a[0] - min_a[0])
      a_d[d][1] = (a_d[d][1] - min_a[1]) / (max_a[1] - min_a[1])
      a_d[d][2] = (a_d[d][2] - min_a[2]) / (max_a[2] - min_a[2])

      # CNN feature normalization (CNN特徴の正規化)
      if CNNmode == 1:
        for m in xrange(M[d]):
          for i in xrange(dim_o):
            o_dm[d][m][i] = (o_dm[d][m][i] - min_o) / (max_o - min_o)

    ###test data###
    #D = 10
    #M = [1 for d in xrange(D)]
    #N = [4 for d in xrange(D)]
    #Ad = [0 for d in xrange(D)]
    #basyo = [[],[-10+uniform(-1,1),-10+uniform(-1,1)], [10+uniform(-1,1),10+uniform(-1,1)], [0+uniform(-1,1),0+uniform(-1,1)]]
    #a_d  = [ basyo[1] , basyo[2] , basyo[3] , basyo[1] , basyo[2] , basyo[3] , basyo[1] , basyo[2] ]#, basyo[3] , basyo[1]]
    #p_dm = [[basyo[3]],[basyo[3]],[basyo[3]],[basyo[1]],[basyo[1]],[basyo[1]],[basyo[2]],[basyo[2]],[basyo[2]],[basyo[3]]]
    #o_dm = [[basyo[1]],[basyo[3]],[basyo[2]],[basyo[3]],[basyo[2]],[basyo[1]],[basyo[2]],[basyo[3]],[basyo[1]],[basyo[2]]]
    #c_dm = [[basyo[2]],[basyo[1]],[basyo[1]],[basyo[2]],[basyo[3]],[basyo[2]],[basyo[3]],[basyo[1]],[basyo[3]],[basyo[1]]]

    #w_dn = [["reach","front","box","green"] for i in range(1)]+[["touch","right","cup","green"] for i in range(1)]+[["touch","front","box","red"] for i in range(1)]+[["reach","front","box","blue"] for i in range(1)]+[["touch","front","box","green"] for i in range(1)]+[["lookat","front","cup","blue"] for i in range(1)]+[["lookat","left","cup","red"] for i in range(1)]+[["lookat","right","box","red"] for i in range(1)]#+[["a3a","p2p","o1o","c3c"] for i in range(1)]+[["a1a","p3p","o2o","c1c"] for i in range(1)]

    # Save the loaded data to one file (読み込んだデータを保存)
    fp = open( foldername +"("+str(sn).zfill(3)+"-"+str(en).zfill(3)+")"+'/' + filename +'/'+ trialname + '_' + filename +'_data.csv', 'w')
    for i in xrange(3):
      fp.write(repr(min_a[i])+',')
    fp.write('\n')
    for i in xrange(3):
      fp.write(repr(max_a[i])+',')
    fp.write('\n')
    fp.write('M\n')
    fp.write(repr(M))
    fp.write('\n')
    fp.write('N\n')
    fp.write(repr(N))
    fp.write('\n')
    fp.write('a_d\n')
    for d in xrange(D):
      fp.write(repr(a_d[d])+'\n')
    fp.write('\n')
    fp.write('p_dm\n')
    fp.write(repr(p_dm))
    fp.write('\n')
    fp.write('o_dm\n')
    for d in xrange(D):
      fp.write(repr(o_dm[d])+'\n')
    fp.write('\n')
    fp.write('c_dm\n')
    fp.write(repr(c_dm))
    fp.write('\n')
    fp.write('Ad\n')
    fp.write(repr(Ad))
    fp.write('\n')
    fp.write('w_dn\n')
    fp.write(repr(w_dn))
    fp.write('\n')
    fp.close()

    return M, N, w_dn, a_d, p_dm, o_dm, c_dm, Ad

def para_save(foldername,trialname,filename,za,zp,zo,zc,Fd,theta,W_list,Mu_a,Sig_a,Mu_p,Sig_p,Mu_o,Sig_o,Mu_c,Sig_c,pi_a,pi_p,pi_o,pi_c):
    foldername = foldername + "/" + filename
    trialname = trialname + "_" + filename

    #Saving each parameter
    fp = open( foldername +'/' + trialname +'_kekka.csv', 'w')
    fp.write('sampling_data\n')
    fp.write('za\n')
    for d in xrange(D):
      fp.write(repr(d)+',')
    fp.write('\n')
    for d in xrange(D):
      fp.write(repr(za[d])+',')
    fp.write('\n')

    fp.write('zp\n')
    for d in xrange(D):
      for m in xrange(len(zp[d])):
        fp.write(repr(d)+'->'+repr(m)+',')
    fp.write('\n')
    for d in xrange(D):
      for m in xrange(len(zp[d])):
        fp.write(repr(zp[d][m])+',')
    fp.write('\n')

    fp.write('zo\n')
    for d in xrange(D):
      for m in xrange(len(zo[d])):
        fp.write(repr(d)+'->'+repr(m)+',')
    fp.write('\n')
    for d in xrange(D):
      for m in xrange(len(zo[d])):
        fp.write(repr(zo[d][m])+',')
    fp.write('\n')

    fp.write('zc\n')
    for d in xrange(D):
      for m in xrange(len(zc[d])):
        fp.write(repr(d)+'->'+repr(m)+',')
    fp.write('\n')
    for d in xrange(D):
      for m in xrange(len(zc[d])):
        fp.write(repr(zc[d][m])+',')
    fp.write('\n')

    fp.write('Fd\n')
    for d in xrange(D):
      fp.write(repr(d)+','+repr(Fd[d])+'\n')

    fp.write('theta\n,,')
    for w in xrange(len(W_list)):
      fp.write(repr(W_list[w])+',')
    fp.write('\n')
    for i in xrange(Ka):
      fp.write('a '+repr(i)+','+repr(i)+',')
      for w in xrange(len(W_list)):
        fp.write(repr(theta[i][w])+',')
      fp.write('\n')
    for i in xrange(Kp):
      fp.write('p '+repr(i)+','+repr(i+dict["p"])+',')
      for w in xrange(len(W_list)):
        fp.write(repr(theta[i+dict["p"]][w])+',')
      fp.write('\n')
    for i in xrange(Ko):
      fp.write('o '+repr(i)+','+repr(i+dict["o"])+',')
      for w in xrange(len(W_list)):
        fp.write(repr(theta[i+dict["o"]][w])+',')
      fp.write('\n')
    for i in xrange(Kc):
      fp.write('c '+repr(i)+','+repr(i+dict["c"])+',')
      for w in xrange(len(W_list)):
        fp.write(repr(theta[i+dict["c"]][w])+',')
      fp.write('\n')

    fp.write('action category\n')
    fp.write('Mu\n')
    for k in xrange(Ka):
      fp.write(repr(k)+',')
      for dim in xrange(dim_a):
        fp.write(repr(Mu_a[k][dim])+',')
      fp.write('\n')
    fp.write('Sig\n')
    for k in xrange(Ka):
      fp.write(repr(k)+'\n')
      for dim in xrange(dim_a):
        for dim2 in xrange(dim_a):
          fp.write(repr(Sig_a[k][dim][dim2])+',')
        fp.write('\n')
      fp.write('\n')

    fp.write('position category\n')
    fp.write('Mu\n')
    for k in xrange(Kp):
      fp.write(repr(k)+',')
      for dim in xrange(dim_p):
        fp.write(repr(Mu_p[k][dim])+',')
      fp.write('\n')
    fp.write('Sig\n')
    for k in xrange(Kp):
      fp.write(repr(k)+'\n')
      for dim in xrange(dim_p):
        for dim2 in xrange(dim_p):
          fp.write(repr(Sig_p[k][dim][dim2])+',')
        fp.write('\n')
      fp.write('\n')

    fp.write('object category\n')
    fp.write('Mu\n')
    for k in xrange(Ko):
      fp.write(repr(k)+',')
      for dim in xrange(dim_o):
        fp.write(repr(Mu_o[k][dim])+',')
      fp.write('\n')
    fp.write('Sig\n')
    for k in xrange(Ko):
      fp.write(repr(k)+'\n')
      for dim in xrange(dim_o):
        for dim2 in xrange(dim_o):
          fp.write(repr(Sig_o[k][dim][dim2])+',')
        fp.write('\n')
      fp.write('\n')

    fp.write('color category\n')
    fp.write('Mu\n')
    for k in xrange(Kc):
      fp.write(repr(k)+',')
      for dim in xrange(dim_c):
        fp.write(repr(Mu_c[k][dim])+',')
      fp.write('\n')
    fp.write('Sig\n')
    for k in xrange(Kc):
      fp.write(repr(k)+'\n')
      for dim in xrange(dim_c):
        for dim2 in xrange(dim_c):
          fp.write(repr(Sig_c[k][dim][dim2])+',')
        fp.write('\n')
      fp.write('\n')

    fp.write('pi_a'+',')
    for k in xrange(Ka):
      fp.write(repr(pi_a[k])+',')
    fp.write('\n')
    fp.write('pi_p'+',')
    for k in xrange(Kp):
      fp.write(repr(pi_p[k])+',')
    fp.write('\n')
    fp.write('pi_o'+',')
    for k in xrange(Ko):
      fp.write(repr(pi_o[k])+',')
    fp.write('\n')
    fp.write('pi_c'+',')
    for k in xrange(Kc):
      fp.write(repr(pi_c[k])+',')
    fp.write('\n')

    fp.close()

    #print 'File Output Successful!(filename:'+filename+')\n'

    ## Output each parameter as each file (パラメータそれぞれをそれぞれのファイルとして出力)
    fp = open(foldername +'/' + trialname +'_za.csv', 'w')
    for d in xrange(D):
      fp.write(repr(za[d])+',')
    fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_zp.csv', 'w')
    for d in xrange(D):
      for m in xrange(len(zp[d])):
        fp.write(repr(zp[d][m])+',')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_zo.csv', 'w')
    for d in xrange(D):
      for m in xrange(len(zo[d])):
        fp.write(repr(zo[d][m])+',')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_zc.csv', 'w')
    for d in xrange(D):
      for m in xrange(len(zc[d])):
        fp.write(repr(zc[d][m])+',')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_Fd.csv', 'w')
    #fp.write('Fd\n')
    for d in xrange(D):
      for f in xrange(len(Fd[d])):
        fp.write(repr(Fd[d][f])+',')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_theta.csv', 'w')
    #fp.write('theta\n')
    for i in xrange(L):
      for w in xrange(len(W_list)):
        fp.write(repr(theta[i][w])+',')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_W_list.csv', 'w')
    for w in xrange(len(W_list)):
      fp.write(repr(W_list[w])+',')
    fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_Mu_a.csv', 'w')
    #fp.write('action category\n')
    #fp.write('Mu\n')
    for k in xrange(Ka):
      #fp.write(repr(k)+',')
      for dim in xrange(dim_a):
        fp.write(repr(Mu_a[k][dim])+',')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_Sig_a.csv', 'w')
    #fp.write('Sig\n')
    for k in xrange(Ka):
      #fp.write(repr(k)+',')
      for dim in xrange(dim_a):
        for dim2 in xrange(dim_a):
          fp.write(repr(Sig_a[k][dim][dim2])+',')
        fp.write('\n')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_Mu_p.csv', 'w')
    #fp.write('position category\n')
    #fp.write('Mu\n')
    for k in xrange(Kp):
      #fp.write(repr(k)+',')
      for dim in xrange(dim_p):
        fp.write(repr(Mu_p[k][dim])+',')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_Sig_p.csv', 'w')
    #fp.write('Sig\n')
    for k in xrange(Kp):
      #fp.write(repr(k)+',')
      for dim in xrange(dim_p):
        for dim2 in xrange(dim_p):
          fp.write(repr(Sig_p[k][dim][dim2])+',')
        fp.write('\n')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_Mu_o.csv', 'w')
    #fp.write('object category\n')
    #fp.write('Mu\n')
    for k in xrange(Ko):
      #fp.write(repr(k)+',')
      for dim in xrange(dim_o):
        fp.write(repr(Mu_o[k][dim])+',')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_Sig_o.csv', 'w')
    #fp.write('Sig\n')
    for k in xrange(Ko):
      #fp.write(repr(k)+',')
      for dim in xrange(dim_o):
        for dim2 in xrange(dim_o):
          fp.write(repr(Sig_o[k][dim][dim2])+',')
        fp.write('\n')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_Mu_c.csv', 'w')
    #fp.write('color category\n')
    #fp.write('Mu\n')
    for k in xrange(Kc):
      #fp.write(repr(k)+',')
      for dim in xrange(dim_c):
        fp.write(repr(Mu_c[k][dim])+',')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_Sig_c.csv', 'w')
    #fp.write('Sig\n')
    for k in xrange(Kc):
      #fp.write(repr(k)+',')
      for dim in xrange(dim_c):
        for dim2 in xrange(dim_c):
          fp.write(repr(Sig_c[k][dim][dim2])+',')
        fp.write('\n')
      fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_pi_a.csv', 'w')
    #fp.write('pi_a'+',')
    for k in xrange(Ka):
      fp.write(repr(pi_a[k])+',')
    fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_pi_p.csv', 'w')
    #fp.write('pi_p'+',')
    for k in xrange(Kp):
      fp.write(repr(pi_p[k])+',')
    fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_pi_o.csv', 'w')
    #fp.write('pi_o'+',')
    for k in xrange(Ko):
      fp.write(repr(pi_o[k])+',')
    fp.write('\n')
    fp.close()

    fp = open(foldername +'/' + trialname +'_pi_c.csv', 'w')
    #fp.write('pi_c'+',')
    for k in xrange(Kc):
      fp.write(repr(pi_c[k])+',')
    fp.write('\n')
    fp.close()

    print 'File Output Successful!(filename:'+foldername+')\n'

# Simulation
def simulate(foldername,trialname,filename,sn,en, M, N, w_dn, a_d, p_dm, o_dm, c_dm, Ad):
      np.random.seed()
      print w_dn
      # random initialzation
      print u"Initialize Parameters..."
      za = [ int(uniform(0,Ko)) for d in xrange(D) ]
      zp = [ [ int(uniform(0,Kp)) for m in xrange(M[d]) ] for d in xrange(D) ]
      zo = [ [ int(uniform(0,Ko)) for m in xrange(M[d]) ] for d in xrange(D) ]
      zc = [ [ int(uniform(0,Kc)) for m in xrange(M[d]) ] for d in xrange(D) ]

      Fd = [ Sample_Frame(N[d]) for d in xrange(D)] #random initialzation ("a","p","o","c")

      cw = np.sum([collections.Counter(w_dn[d]) for d in xrange(D)])
      W_list = list(cw)  ##List of words
      W = len(cw)  ##Length of word list
      theta = [ sum(dirichlet(np.array([gamma for w in xrange(W)]),100))/100.0 for i in xrange(L) ] #indexと各モダリティーの対応付けはdictionary形式で呼び出す

      #KMeans(n_clusters=Ka, init='k-means++').fit(a_d).cluster_centers_
      Mu_a  = KMeans(n_clusters=Ka, init='k-means++').fit(a_d).cluster_centers_#[ np.array([uniform(mu_a_init[0],mu_a_init[1]) for i in xrange(dim_a)]) for k in xrange(Ka) ]
      Sig_a = [ np.eye(dim_a)*sig_a_init for k in xrange(Ka) ]
      #print "Mu_a",Mu_a

      p_temp = []
      for d in xrange(D):
        p_temp = p_temp + p_dm[d]
      Mu_p  = KMeans(n_clusters=Kp, init='k-means++').fit(p_temp).cluster_centers_#[ np.array([uniform(mu_p_init[0],mu_p_init[1]) for i in xrange(dim_p)]) for k in xrange(Kp) ]
      Sig_p = [ np.eye(dim_p)*sig_p_init for k in xrange(Kp) ]

      o_temp = []
      for d in xrange(D):
        o_temp = o_temp + o_dm[d]
      Mu_o  = KMeans(n_clusters=Ko, init='k-means++').fit(o_temp).cluster_centers_#[ np.array([uniform(mu_o_init[0],mu_o_init[1]) for i in xrange(dim_o)]) for k in xrange(Ko) ]
      Sig_o = [ np.eye(dim_o)*sig_o_init for k in xrange(Ko) ]

      c_temp = []
      for d in xrange(D):
        c_temp = c_temp + c_dm[d]
      Mu_c  = KMeans(n_clusters=Kc, init='k-means++').fit(c_temp).cluster_centers_#[ np.array([uniform(mu_c_init[0],mu_c_init[1]) for i in xrange(dim_c)]) for k in xrange(Kc) ]
      Sig_c = [ np.eye(dim_c)*sig_c_init for k in xrange(Kc) ]

      if nonpara == 0 :
        pi_a = sum(dirichlet([ alpha_a for c in xrange(Ka)],100))/100.0
        pi_p = sum(dirichlet([ alpha_p for c in xrange(Kp)],100))/100.0
        pi_o = sum(dirichlet([ alpha_o for c in xrange(Ko)],100))/100.0
        pi_c = sum(dirichlet([ alpha_c for c in xrange(Kc)],100))/100.0
      elif nonpara == 1:
        pi_a = stick_breaking(alpha_a, Ka)
        pi_p = stick_breaking(alpha_p, Kp)
        pi_o = stick_breaking(alpha_o, Ko)
        pi_c = stick_breaking(alpha_c, Kc)

      print theta
      print pi_a
      print pi_p
      print pi_o
      print pi_c
      print Mu_a
      print Mu_p
      print Mu_o
      print Mu_c
      print Fd

      ### Copy initial values (このやり方でないと値が変わってしまう)
      za_init = [ za[d] for d in xrange(D)]
      zp_init = [ [ zp[d][m] for m in xrange(M[d]) ] for d in xrange(D) ]
      zo_init = [ [ zo[d][m] for m in xrange(M[d]) ] for d in xrange(D) ]
      zc_init = [ [ zc[d][m] for m in xrange(M[d]) ] for d in xrange(D) ]
      Fd_init = [ Fd[d] for d in xrange(D)]
      theta_init = [ np.array(theta[i]) for i in xrange(L) ]
      Mu_a_init  = [ np.array(Mu_a[k]) for k in xrange(Ka) ]
      Sig_a_init = [ np.array(Sig_a[k]) for k in xrange(Ka) ]
      Mu_p_init  = [ np.array(Mu_p[k]) for k in xrange(Kp) ]
      Sig_p_init = [ np.array(Sig_p[k]) for k in xrange(Kp) ]
      Mu_o_init  = [ np.array(Mu_o[k]) for k in xrange(Ko) ]
      Sig_o_init = [ np.array(Sig_o[k]) for k in xrange(Ko) ]
      Mu_c_init  = [ np.array(Mu_c[k]) for k in xrange(Kc) ]
      Sig_c_init = [ np.array(Sig_c[k]) for k in xrange(Kc) ]
      pi_a_init = [pi_a[k] for k in xrange(Ka)]
      pi_p_init = [pi_p[k] for k in xrange(Kp)]
      pi_o_init = [pi_o[k] for k in xrange(Ko)]
      pi_c_init = [pi_c[k] for k in xrange(Kc)]

      ### save initial values
      #filename_init = filename + "/init"
      trialname_init = "init/" + trialname
      para_save(foldername,trialname_init,filename,za_init,zp_init,zo_init,zc_init,Fd_init,theta_init,W_list,Mu_a_init,Sig_a_init,Mu_p_init,Sig_p_init,Mu_o_init,Sig_o_init,Mu_c_init,Sig_c_init,pi_a_init,pi_p_init,pi_o_init,pi_c_init)

      ######################################################################
      ####                     ↓ Learning phase ↓                       ####
      ######################################################################
      print u"- <START> Learning of Lexicon and Multiple Categories ver. iCub MODEL. -"

      for iter in xrange(num_iter):   #Iteration of Gibbs sampling
        print '----- Iter. '+repr(iter+1)+' -----'

        ########## ↓ ##### Sampling za ##### ↓ ##########
        print u"Sampling za..."

        for d in xrange(D):         #for each sentence
          temp = np.array(pi_a)
          for k in xrange(Ka):      #for each index of category
            for n in xrange(N[d]):  #for each word in a sentence
              if Fd[d][n] == "a":
                temp[k] = temp[k] * theta[k + dict["a"]][W_list.index(w_dn[d][n])]
            temp[k] = temp[k] * multivariate_normal.pdf(a_d[d], mean=Mu_a[k], cov=Sig_a[k])

          temp = temp / np.sum(temp)  #正規化
          za[d] = list(multinomial(1,temp)).index(1)
        print za
        ########## ↑ ##### Sampling za ##### ↑ ##########

        ########## ↓ ##### Sampling zp ##### ↓ ##########
        print u"Sampling zp..."

        for d in xrange(D):         #for each sentence
          for m in xrange(M[d]):    #for each obect
            temp = np.array(pi_p)
            for k in xrange(Kp):      #for each index of category
              for n in xrange(N[d]):  #for each word in a sentence
                if Fd[d][n] == "p" and Ad[d] == m:
                  temp[k] = temp[k] * theta[k + dict["p"]][W_list.index(w_dn[d][n])]
              temp[k] = temp[k] * multivariate_normal.pdf(p_dm[d][m], mean=Mu_p[k], cov=Sig_p[k])

            temp = temp / np.sum(temp)  #正規化
            zp[d][m] = list(multinomial(1,temp)).index(1)
        print zp
        ########## ↑ ##### Sampling zp ##### ↑ ##########

        ########## ↓ ##### Sampling zo ##### ↓ ##########
        print u"Sampling zo..."

        for d in xrange(D):         #データごとに
          for m in xrange(M[d]):    #物体ごとに
            temp = np.array(pi_o)
            logtemp = np.array([log(pi_o[k]) for k in xrange(Ko)])
            for k in xrange(Ko):      #カテゴリ番号ごとに
              for n in xrange(N[d]):  #文中の単語ごとに
                if Fd[d][n] == "o" and Ad[d] == m:
                  #temp[k] = temp[k] * theta[k + dict["o"]][W_list.index(w_dn[d][n])]
                  logtemp[k] = logtemp[k] + log(theta[k + dict["o"]][W_list.index(w_dn[d][n])])
              loggauss = multivariate_normal.logpdf(o_dm[d][m], mean=Mu_o[k], cov=Sig_o[k])
              #print loggauss
              logtemp[k] = logtemp[k] + loggauss#multivariate_normal.pdf(o_dm[d][m], mean=Mu_o[k], cov=Sig_o[k])

            logtemp = logtemp - np.max(logtemp)
            logtemp = logtemp - sp.misc.logsumexp(logtemp)#log(np.sum(np.array([exp(logtemp[k]) for k in xrange(Ko)])))  #正規化
            #print logtemp,sp.misc.logsumexp(logtemp)
            #print np.exp(logtemp),np.sum(np.exp(logtemp))
            zo[d][m] = list( multinomial(1,np.exp(logtemp)) ).index(1)
        print zo
        ########## ↑ ##### Sampling zo ##### ↑ ##########

        ########## ↓ ##### Sampling zc ##### ↓ ##########
        print u"Sampling zc..."

        for d in xrange(D):         #データごとに
          for m in xrange(M[d]):    #物体ごとに
            temp = np.array(pi_c)
            for k in xrange(Kc):      #カテゴリ番号ごとに
              for n in xrange(N[d]):  #文中の単語ごとに
                if Fd[d][n] == "c" and Ad[d] == m:
                  #print temp
                  temp[k] = temp[k] * theta[k + dict["c"]][W_list.index(w_dn[d][n])]
                  #print temp
              temp[k] = temp[k] * multivariate_normal.pdf(c_dm[d][m], mean=Mu_c[k], cov=Sig_c[k])
              #print temp[k],multivariate_normal.pdf(c_dm[d][m], mean=Mu_c[k], cov=Sig_c[k])
            #print temp
            temp = temp / np.sum(temp)  #正規化
            #print temp
            zc[d][m] = list(multinomial(1,temp)).index(1)
        print zc
        ########## ↑ ##### Sampling zc ##### ↑ ##########

        ########## ↓ ##### π_aのサンプリング ##### ↓ ##########
        print u"Sampling PI_a..."

        cc = collections.Counter(za)
        if nonpara == 0:
          temp = np.array([cc[k] + alpha_a for k in xrange(Ka)])
        elif nonpara == 1:
          temp = np.array([cc[k] + (alpha_a / float(Ka)) for k in xrange(Ka)])

        #加算したデータとパラメータから事後分布を計算しサンプリング
        pi_a = dirichlet(temp)
        print pi_a
        ########## ↑ ##### π_aのサンプリング ##### ↑ ##########

        ########## ↓ ##### π_pのサンプリング ##### ↓ ##########
        print u"Sampling PI_p..."

        cc = np.sum([collections.Counter(zp[d]) for d in range(D)])
        if nonpara == 0:
          temp = np.array([cc[k] + alpha_p for k in xrange(Kp)])
        elif nonpara == 1:
          temp = np.array([cc[k] + (alpha_p / float(Kp)) for k in xrange(Ka)])

        #加算したデータとパラメータから事後分布を計算しサンプリング
        pi_p = dirichlet(temp)
        print pi_p
        ########## ↑ ##### π_pのサンプリング ##### ↑ ##########

        ########## ↓ ##### π_oのサンプリング ##### ↓ ##########
        print u"Sampling PI_o..."

        cc = np.sum([collections.Counter(zo[d]) for d in range(D)])
        if nonpara == 0:
          temp = np.array([cc[k] + alpha_o for k in xrange(Ko)])
        elif nonpara == 1:
          temp = np.array([cc[k] + (alpha_o / float(Ko)) for k in xrange(Ka)])

        #加算したデータとパラメータから事後分布を計算しサンプリング
        pi_o = dirichlet(temp)
        print pi_o
        ########## ↑ ##### π_oのサンプリング ##### ↑ ##########

        ########## ↓ ##### π_cのサンプリング ##### ↓ ##########
        print u"Sampling PI_c..."

        cc = np.sum([collections.Counter(zc[d]) for d in range(D)])
        if nonpara == 0:
          temp = np.array([cc[k] + alpha_c for k in xrange(Kc)])
        elif nonpara == 1:
          temp = np.array([cc[k] + (alpha_c / float(Kc)) for k in xrange(Ka)])

        #加算したデータとパラメータから事後分布を計算しサンプリング
        pi_c = dirichlet(temp)
        print pi_c
        ########## ↑ ##### π_cのサンプリング ##### ↑ ##########

        ########## ↓ ##### μa,Σa(ガウス分布の平均、共分散行列)のサンプリング ##### ↓ ##########
        print u"Sampling myu_a,Sigma_a..."

        cc = collections.Counter(za)
        for k in xrange(Ka) :
          nk = cc[k]
          xt = []
          m_ML = np.zeros(dim_a)
          ###kについて、zaが同じものを集める
          if nk != 0 :  #もしzaの中にkがあれば(計算短縮処理)        ##0ワリ回避
            for d in xrange(D) :
              if za[d] == k :
                xt = xt + [ np.array(a_d[d]) ]

            m_ML = sum(xt) / float(nk) #fsumではダメ
            #print "n:%d m_ML:%s" % (nk,str(m_ML))
            print "a%d n:%d" % (k,nk)

            ##ハイパーパラメータ更新
            kN = k0a + nk
            mN = ( k0a*m0a + nk*m_ML ) / kN  #dim_a 次元横ベクトル
            nN = n0a + nk
            VN = V0a + sum([np.dot(np.array([xt[j]-m_ML]).T,np.array([xt[j]-m_ML])) for j in xrange(nk)]) + (k0a*nk/kN)*np.dot(np.array([m_ML-m0a]).T,np.array([m_ML-m0a]))

            ##3.1##Σを逆ウィシャートからサンプリング
            Sig_a[k] = invwishart.rvs(df=nN, scale=VN) #/ n0a
            ##3.2##μを多変量ガウスからサンプリング
            Mu_a[k] = np.mean([multivariate_normal.rvs(mean=mN, cov=Sig_a[k]/kN) for i in xrange(100)],0) #サンプリングをロバストに
          else:  #データがないとき
            Mu_a[k]  = np.array([uniform(mu_a_init[0],mu_a_init[1]) for i in xrange(dim_a)])
            Sig_a[k] = invwishart.rvs(df=n0a, scale=V0a )#np.eye(dim_a)*sig_a_init

          if (nk != 0):  #データなしは表示しない
            print 'Mu_a '+str(k)+' : '+str(Mu_a[k])
            print 'Sig_a'+str(k)+':\n'+str(Sig_a[k])

          #print [(np.array([xt[j]-m_ML]).T,np.array([xt[j]-m_ML])) for j in xrange(nk)]
          #print [np.dot(np.array([xt[j]-m_ML]).T,np.array([xt[j]-m_ML])) for j in xrange(nk)]
          #print sum([np.dot(np.array([xt[j]-m_ML]).T,np.array([xt[j]-m_ML])) for j in xrange(nk)])
          #print ( float(k0a*nk)/kN ) * np.dot(np.array([m_ML - m0a]).T,np.array([m_ML - m0a]))
          #print VN
          #samp_sig_rand = np.array([ invwishart(nuN,VN) for i in xrange(100)])    ######
          #samp_sig = np.mean(samp_sig_rand,0)
          #print samp_sig

          #if np.linalg.det(samp_sig) < -0.0:
          #  Sig_a[k] = np.mean(np.array([ invwishartrand(nuN,VN)]),0)

          #print ''
          #for k in xrange(K):
          #if (nk[k] != 0):  #データなしは表示しない
          #  print 'Sig_a'+str(k)+':'+str(Sig_a[k])

        ########## ↑ ##### μa,Σa(ガウス分布の平均、共分散行列)のサンプリング ##### ↑ ##########

        ########## ↓ ##### μp,Σp(ガウス分布の平均、共分散行列)のサンプリング ##### ↓ ##########
        print u"Sampling myu_p,Sigma_p..."

        cc = np.sum([collections.Counter(zp[d]) for d in range(D)])
        for k in xrange(Kp) :
          nk = cc[k]
          xt = []
          m_ML = np.zeros(dim_p)
          ###kについて、zaが同じものを集める
          if nk != 0 :  #もしzaの中にkがあれば(計算短縮処理)        ##0ワリ回避
            for d in xrange(D) :
              for m in xrange(M[d]):
                if zp[d][m] == k :
                  xt = xt + [ np.array(p_dm[d][m]) ]

            m_ML = sum(xt) / float(nk) #fsumではダメ
            #print "n:%d m_ML:%s" % (nk,str(m_ML))
            print "p%d n:%d" % (k,nk)

            ##ハイパーパラメータ更新
            kN = k0p + nk
            mN = ( k0p*m0p + nk*m_ML ) / kN  #dim_a 次元横ベクトル
            nN = n0p + nk
            VN = V0p + sum([np.dot(np.array([xt[j]-m_ML]).T,np.array([xt[j]-m_ML])) for j in xrange(nk)]) + (k0p*nk/kN)*np.dot(np.array([m_ML-m0p]).T,np.array([m_ML-m0p]))

            ##3.1##Σを逆ウィシャートからサンプリング
            Sig_p[k] = invwishart.rvs(df=nN, scale=VN) #/ n0a
            ##3.2##μを多変量ガウスからサンプリング
            Mu_p[k] = np.mean([multivariate_normal.rvs(mean=mN, cov=Sig_p[k]/kN) for i in xrange(100)],0) #サンプリングをロバストに
          else:  #データがないとき
            Mu_p[k]  = np.array([uniform(mu_p_init[0],mu_p_init[1]) for i in xrange(dim_p)])
            Sig_p[k] = invwishart.rvs(df=n0p, scale=V0p ) #np.eye(dim_p)*sig_p_init

          if (nk != 0):  #データなしは表示しない
            print 'Mu_p '+str(k)+' : '+str(Mu_p[k])
            print 'Sig_p'+str(k)+':\n'+str(Sig_p[k])

        ########## ↑ ##### μp,Σp(ガウス分布の平均、共分散行列)のサンプリング ##### ↑ ##########

        ########## ↓ ##### μo,Σo(ガウス分布の平均、共分散行列)のサンプリング ##### ↓ ##########
        print u"Sampling myu_o,Sigma_o..."

        cc = np.sum([collections.Counter(zo[d]) for d in range(D)])
        for k in xrange(Ko) :
          nk = cc[k]
          xt = []
          m_ML = np.zeros(dim_o)
          ###kについて、zaが同じものを集める
          if nk != 0 :  #もしzaの中にkがあれば(計算短縮処理)        ##0ワリ回避
            for d in xrange(D) :
              for m in xrange(M[d]):
                if zo[d][m] == k :
                  xt = xt + [ np.array(o_dm[d][m]) ]

            m_ML = sum(xt) / float(nk) #fsumではダメ
            #print "n:%d m_ML:%s" % (nk,str(m_ML))
            print "o%d n:%d" % (k,nk)

            ##ハイパーパラメータ更新
            kN = k0o + nk
            mN = ( k0o*m0o + nk*m_ML ) / kN  #dim_a 次元横ベクトル
            nN = n0o + nk
            VN = V0o + sum([np.dot(np.array([xt[j]-m_ML]).T,np.array([xt[j]-m_ML])) for j in xrange(nk)]) + (k0o*nk/kN)*np.dot(np.array([m_ML-m0o]).T,np.array([m_ML-m0o]))

            ##3.1##Σを逆ウィシャートからサンプリング
            Sig_o[k] = invwishart.rvs(df=nN, scale=VN) #/ n0a
            ##3.2##μを多変量ガウスからサンプリング
            Mu_o[k] = np.mean([multivariate_normal.rvs(mean=mN, cov=Sig_o[k]/kN) for i in xrange(100)],0) #サンプリングをロバストに
          else:  #データがないとき
            Mu_o[k]  = np.array([uniform(mu_o_init[0],mu_o_init[1]) for i in xrange(dim_o)])
            Sig_o[k] = invwishart.rvs(df=n0o, scale=V0o ) #np.eye(dim_o)*sig_o_init

          if (nk != 0):  #データなしは表示しない
            print 'Mu_o '+str(k)+' : '+str(Mu_o[k])
            print 'Sig_o'+str(k)+':\n'+str(Sig_o[k])

        ########## ↑ ##### μo,Σo(ガウス分布の平均、共分散行列)のサンプリング ##### ↑ ##########

        ########## ↓ ##### μc,Σc(ガウス分布の平均、共分散行列)のサンプリング ##### ↓ ##########
        print u"Sampling myu_c,Sigma_c..."

        cc = np.sum([collections.Counter(zc[d]) for d in range(D)])
        for k in xrange(Kc) :
          nk = cc[k]
          xt = []
          m_ML = np.zeros(dim_c)
          ###kについて、zaが同じものを集める
          if nk != 0 :  #もしzaの中にkがあれば(計算短縮処理)        ##0ワリ回避
            for d in xrange(D) :
              for m in xrange(M[d]):
                if zc[d][m] == k :
                  xt = xt + [ np.array(c_dm[d][m]) ]

            m_ML = sum(xt) / float(nk) #fsumではダメ
            #print "n:%d m_ML:%s" % (nk,str(m_ML))
            print "c%d n:%d" % (k,nk)

            ##ハイパーパラメータ更新
            kN = k0c + nk
            mN = ( k0c*m0c + nk*m_ML ) / kN  #dim_a 次元横ベクトル
            nN = n0c + nk
            VN = V0c + sum([np.dot(np.array([xt[j]-m_ML]).T,np.array([xt[j]-m_ML])) for j in xrange(nk)]) + (k0c*nk/kN)*np.dot(np.array([m_ML-m0c]).T,np.array([m_ML-m0c]))

            ##3.1##Σを逆ウィシャートからサンプリング
            Sig_c[k] = invwishart.rvs(df=nN, scale=VN) #/ n0a
            ##3.2##μを多変量ガウスからサンプリング
            Mu_c[k] = np.mean([multivariate_normal.rvs(mean=mN, cov=Sig_c[k]/kN) for i in xrange(100)],0) #サンプリングをロバストに
          else:  #データがないとき
            Mu_c[k]  = np.array([uniform(mu_c_init[0],mu_c_init[1]) for i in xrange(dim_c)])
            Sig_c[k] = invwishart.rvs(df=n0c, scale=V0c ) #np.eye(dim_c)*sig_c_init

          if (nk != 0):  #データなしは表示しない
            print 'Mu_c '+str(k)+' : '+str(Mu_c[k])
            print 'Sig_c'+str(k)+':\n'+str(Sig_c[k])

        ########## ↑ ##### μc,Σc(ガウス分布の平均、共分散行列)のサンプリング ##### ↑ ##########

        ########## ↓ ##### Fdのサンプリング ##### ↓ ##########
        print u"Sampling Fd..."
        #基本的にデータごとの単語数におけるモダリティの組み合わせ全探索する

        for d in xrange(D):
          F_temp = [f for f in itertools.permutations(modality,N[d])]  ##モダリティの順列組み合わせ
          temp = [1.0 for i in xrange(len(F_temp))]  ##フレームの組み合わせ数分の要素を用意する
          for i in xrange(len(F_temp)):
            for n in xrange(N[d]):
              #print i,n,N[d],M[d]
              if F_temp[i][n] == "a":
                temp[i] = temp[i] * theta[za[d]                   ][W_list.index(w_dn[d][n])]
              if F_temp[i][n] == "p":
                temp[i] = temp[i] * theta[zp[d][Ad[d]] + dict["p"]][W_list.index(w_dn[d][n])]
              if F_temp[i][n] == "o":
                temp[i] = temp[i] * theta[zo[d][Ad[d]] + dict["o"]][W_list.index(w_dn[d][n])]
              if F_temp[i][n] == "c":
                temp[i] = temp[i] * theta[zc[d][Ad[d]] + dict["c"]][W_list.index(w_dn[d][n])]
          temp = temp / np.sum(temp)  #正規化
          Fd[d] = F_temp[list(multinomial(1,temp)).index(1)]
          print d, Fd[d]

        ########## ↑ ##### Fdのサンプリング ##### ↑ ##########

        ########## ↓ ##### Sampling Θ ##### ↓ ##########
        print u"Sampling Theta..."
        #dict = {"a":0, "p":Ka, "o":Ka+Kp, "c":Ka+Kp+Ko}   #各モダリティのindexにキーを足すとΘのindexになる

        temp = [np.array([gamma for w in xrange(W)]) for i in xrange(L)]
        for d in xrange(D):
            for n in xrange(N[d]):
              if (Fd[d][n] == "a"):
                temp[za[d]]                   [W_list.index(w_dn[d][n])] = temp[za[d]]                   [W_list.index(w_dn[d][n])] + 1
              if (Fd[d][n] == "p"):
                temp[zp[d][Ad[d]] + dict["p"]][W_list.index(w_dn[d][n])] = temp[zp[d][Ad[d]] + dict["p"]][W_list.index(w_dn[d][n])] + 1
              if (Fd[d][n] == "o"):
                temp[zo[d][Ad[d]] + dict["o"]][W_list.index(w_dn[d][n])] = temp[zo[d][Ad[d]] + dict["o"]][W_list.index(w_dn[d][n])] + 1
              if (Fd[d][n] == "c"):
                temp[zc[d][Ad[d]] + dict["c"]][W_list.index(w_dn[d][n])] = temp[zc[d][Ad[d]] + dict["c"]][W_list.index(w_dn[d][n])] + 1

        #Sampling from the posterior distribution calculated by the added data and hyperparameter
        theta = [sum(dirichlet(temp[i],100))/100.0 for i in xrange(L)] ##ロバストネスを上げる100

        print theta
        ########## ↑ ##### Sampling Θ ##### ↑ ##########
        print ""

      ######################################################################
      ####                     ↑ Learning phase ↑                       ####
      ######################################################################

      loop = 1
      ########  ↓Files output↓  ########
      if loop == 1:
        print "--------------------"
        #最終学習結果を出力
        print u"- <COMPLETED> Learning of Lexicon and Multiple Categories ver. iCub MODEL. -"
        #print 'Sample: ' + str(sample)
        print 'za: ' + str(za)
        print 'zp: ' + str(zp)
        print 'zo: ' + str(zo)
        print 'zc: ' + str(zc)
        for d in xrange(D):
          print 'Fd%d: %s' % (d, str(Fd[d]))
        for c in xrange(Ka):
          print "theta_a%d: %s" % (c,theta[c + dict["a"]])
        for c in xrange(Kp):
          print "theta_p%d: %s" % (c,theta[c + dict["p"]])
        for c in xrange(Ko):
          print "theta_o%d: %s" % (c,theta[c + dict["o"]])
        for c in xrange(Kc):
          print "theta_c%d: %s" % (c,theta[c + dict["c"]])
        for k in xrange(Ka):
          print "mu_a%d: %s" % (k, str(Mu_a[k]))
        for k in xrange(Kp):
          print "mu_p%d: %s" % (k, str(Mu_p[k]))
        for k in xrange(Ko):
          print "mu_o%d: %s" % (k, str(Mu_o[k]))
        for k in xrange(Kc):
          print "mu_c%d: %s" % (k, str(Mu_c[k]))
        #for k in xrange(K):
        #  print "sig%d: \n%s" % (k, str(Sig_a[k]))
        print 'pi_a: ' + str(pi_a)
        print 'pi_p: ' + str(pi_p)
        print 'pi_o: ' + str(pi_o)
        print 'pi_c: ' + str(pi_c)

        print "--------------------"

        #Saving parameters to files
        para_save(foldername,trialname,filename,za,zp,zo,zc,Fd,theta,W_list,Mu_a,Sig_a,Mu_p,Sig_p,Mu_o,Sig_o,Mu_c,Sig_c,pi_a,pi_p,pi_o,pi_c)

      ########  ↑Files output↑  ########

if __name__ == '__main__':
    import sys
    import shutil
    from __init__ import *

    #filename = sys.argv[1]
    #print filename

    #出力ファイル名を要求
    #filename = raw_input("trialname?(folder) >")

    trialname = raw_input("trialname?(folder) >")
    start = raw_input("start number?>")
    end   = raw_input("end number?>")
    filename = raw_input("learning trial name?>")

    sn = int(start)
    en = int(end)
    Data = int(en) - int(sn) +1

    if D != Data:
      print "data number error.",D,Data
      exit()

    foldername = datafolder + trialname+"("+str(sn).zfill(3)+"-"+str(en).zfill(3)+")"

    #make folder
    Makedir( foldername + "/" + filename )
    Makedir( foldername + "/" + filename + "/init")

    #Copy of init.py
    shutil.copy("./__init__.py", foldername + "/" + filename + "/init")

    #Reading data files (単語、切り出した物体特徴と色と位置、姿勢と触覚と手先位置)
    M, N, w_dn, a_d, p_dm, o_dm, c_dm, Ad = data_read(trialname,filename,sn,en)
    #print w_dn

    #Gibbs sampling
    simulate(foldername,trialname,filename,sn,en, M, N, w_dn, a_d, p_dm, o_dm, c_dm, Ad)

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

max_stars_repo_name: a-taniguchi/CSL-BGM

max_issues_repo_name: a-taniguchi/CSL-BGM

max_forks_repo_name: a-taniguchi/CSL-BGM

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


text:
import numpy as np

def myplus(a, b=0):
    """
    Parameters
    ----------
    a : float
        the first number
    b : float
        the second number
        defaults to zero

    Returns
    -------
    a + b
    """
    return a + b

max_stars_repo_name: hypergravity/csstdemo

max_issues_repo_name: hypergravity/csstdemo

max_forks_repo_name: hypergravity/csstdemo