text:
"""Module for symbolical treatment of smooth reservoir models.

This module handles the symbolic treatment of compartmental/reservoir/pool
models.
It does not deal with numerical computations and model simulations,
but rather defines the underlying structure of the respective model.

All fluxes or matrix entries are supposed to be SymPy expressions.
*Smooth* means that no ``Piecewise`` or ``DiracDelta`` functions should be
involved in the model description.

Counting of compartment/pool/reservoir numbers starts at zero and the
total number of pools is :math:`d`.
"""

import multiprocessing
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import numpy as np
import networkx as nx

from copy import copy, deepcopy
from string import Template
from functools import reduce
from sympy import (zeros, Matrix, simplify, diag, eye, gcd, latex, Symbol,
                   flatten, Function, solve, limit, oo , ask , Q, assuming
                   ,sympify)
from sympy.printing import pprint
from . import helpers_reservoir as hr
from .cs_plotter import CSPlotter
from typing import TypeVar

class Error(Exception):
    """Generic error occurring in this module."""
    pass

class SmoothReservoirModel(object):
    """General class of smooth reservoir models.

    Attributes:
        state_vector (SymPy dx1-matrix): The model's state vector
            :math:`x`.
            Its entries are SymPy symbols.
        state_variables (list of str):
            Names of the variables in the state vector.
            Its entries are of type ``str``.
        time_symbol (SymPy symbol): The model's time symbol.
        input_fluxes (dict): The model's external input fluxes.
            ``{key1: flux1, key2: flux2}`` with ``key`` the pool number and
            ``flux`` a SymPy expression for the influx.
        output_fluxes (dict): The model's external output fluxes.
            ``{key1: flux1, key2: flux2}`` with ``key`` the pool number and
            ``flux`` a SymPy expression for the outflux.
        internal_fluxes (dict): The model's internal_fluxes.
            ``{key1: flux1, key2: flux2}`` with ``key = (pool_from, pool_to)``
            and *flux* a SymPy expression for the flux.
    """
    @classmethod
    def from_state_variable_indexed_fluxes(cls,state_vector, time_symbol,
            input_fluxes, output_fluxes, internal_fluxes)->"SmoothReservoirModel":
        """Return an instance of SmoothReservoirModel.

        Args:
            state_vector (SymPy dx1-matrix): The model's state vector
                :math:`x`.
                Its entries are SymPy symbols.
            time_symbol (SymPy symbol): The model's time symbol.
            input_fluxes (dict): The model's external input fluxes.
                ``{key1: flux1, key2: flux2}`` with ``key`` the symbol of the target pool. (as used in the state vector)
                and ``flux`` a SymPy expression for the influx.
            output_fluxes (dict): The model's external output fluxes.
                ``{key1: flux1, key2: flux2}`` with ``key`` the symbols for the source pool (as used in the state vector)
                and ``flux`` a SymPy expression for the outflux.
            internal_fluxes (dict): The model's internal_fluxes.
                ``{key1: flux1, key2: flux2}`` with
                ``key = (source pool symbol, target pool symbol)`` and ``flux`` a SymPy expression
                for the flux.

        Returns:
            :class:`SmoothReservoirModel`
        """
        # transform to integer indexed dicts
        int_input = hr.to_int_keys_1(input_fluxes, state_vector)
        int_output = hr.to_int_keys_1(output_fluxes, state_vector)
        int_internal = hr.to_int_keys_2(internal_fluxes, state_vector)
        # call normal init
        return cls(state_vector, time_symbol, int_input, int_output, int_internal)

    def __init__(self, state_vector, time_symbol,
                       input_fluxes={}, output_fluxes={}, internal_fluxes={}):
        """Initialize an instance of SmoothReservoirModel.

        Args:
            state_vector (SymPy dx1-matrix): The model's state vector
                :math:`x`.
                Its entries are SymPy symbols.
            time_symbol (SymPy symbol): The model's time symbol.
            input_fluxes (dict): The model's external input fluxes.
                ``{key1: flux1, key2: flux2}`` with ``key`` the pool number
                and ``flux`` a SymPy expression for the influx.
            output_fluxes (dict): The model's external output fluxes.
                ``{key1: flux1, key2: flux2}`` with ``key`` the pool number
                and ``flux`` a SymPy expression for the outflux.
            internal_fluxes (dict): The model's internal_fluxes.
                ``{key1: flux1, key2: flux2}`` with
                ``key = (pool_from, pool_to)`` and ``flux`` a SymPy expression
                for the flux.

        Returns:
            :class:`SmoothReservoirModel`
        """
        self.state_vector = state_vector
        self.state_variables = [sv.name for sv in state_vector]
        self.time_symbol=time_symbol
        self.input_fluxes=input_fluxes
        self.output_fluxes=output_fluxes
        self.internal_fluxes=internal_fluxes

    def is_state_dependent(self,expr):
        efss=expr.free_symbols
        svs=set([e for e in self.state_vector])
        inter=efss.intersection(svs)
        return not(len(inter)==0)

    @property
    def is_linear(self):
        """Returns True if we can make SURE that the model is linear by checking that the jacobian is not state dependent.
        Note that external numerical functions of state variables are represented as sympy.Function f(x_1,x_2,...,t)
        Sympy will consider the derivative of math:`df/dx_i` with respect to state variable math:`x_i` as function math:`g(x_1, x_2,...)` too, since it can not exclude this possibility if we know f only numerically.
        In consequence this method will return False even if the numerical implementation of f IS linear in math:`x_1,x_2,...` .
        To avoid this situation you can just reformulate linear external functions math:`f(x_1,x_2,...,t)` as linear combinations
        of state independent external functions math:`f(x_1,x_2,...,t)=g_1(t)x_1+g_2(t)x_2+...` so that sympy can detect the linearity.

        Returns:
            bool: 'True', 'False'

        """
        return not(self.is_state_dependent(self.jacobian))

    # the following functions are used by the 'figure' method to determine
    # the color of the respective arrow

    # mm 1-17-2022
    # The old implementation using gcd had several issues:
    # 1.) The following code lead to  wrong results
    #     if gcd(sv, flux) == 1:
    #        return 'no state dependence'
    #     gcd(x,sqrt(x))==1 (according to sympy!)
    #     but the flux is obviously nonlinear and state dependent
    # 2.) Furthermore gcd chokes at Piecewise expression
    #     although perfectly sensible (e.g.piecewise in time)

    # a little helper
    def _linear__state_dependent__nonlinear(self, flux, sv):
        # sympify converts constant fluxes (e.g. integers)
        # to sympy expressions that have a free_symbols method
        flux = sympify(flux)

        rate = simplify(flux / sv)
        if all([sv not in rate.free_symbols for sv in list(self.state_vector)]):
            return "linear"
        else:
            # note that no state dependence is a subcase of nonlinear
            if len({sv}.intersection(flux.free_symbols))==0:
                return 'no state dependence'

            else:  #if any([sv in rate.free_symbols for sv in list(state_vector)]):
                return 'nonlinear'

    def _output_flux_type(self, pool_from):
        """Return the type of an external output flux.

        Args:
            pool_from (int): The number of the pool from which the flux starts.

        Returns:
            str: 'linear', 'nonlinear', 'no state dependence'
        """
        sv = self.state_vector[pool_from]

        flux = sympify(self.output_fluxes[pool_from])

        return self._linear__state_dependent__nonlinear(flux, sv)

    def _internal_flux_type(self, pool_from, pool_to):
        """Return the type of an internal flux.

        Args:
            pool_from (int): The number of the pool from which the flux starts.
            pool_to (int): The number of the pool to which the flux goes.

        Returns:
            str: 'linear', 'nonlinear', 'no state dependence'
        """
        sv = self.state_vector[pool_from]
        flux = self.internal_fluxes[(pool_from, pool_to)]
        return self._linear__state_dependent__nonlinear(flux, sv)

    def _input_flux_type(self, pool_to):
        """Return the type of an external input flux.

        Args:
            pool_to (int): The number of the pool to which the flux contributes.

        Returns:
            str: 'linear', 'nonlinear', 'no state dependence'
        """
        sv = self.state_vector[pool_to]
        # we compute the derivative of the appropriate row of the input vector w.r.t. all the state variables
        # (This is a row of the jacobian)
        u_i=Matrix([self.external_inputs[pool_to]])
        s_v=Matrix(self.state_vector)
        J_i=hr.jacobian(u_i,s_v)
        # an input that does not depend on state variables has a zero derivative with respect
        # to all state variables
        if all([ j_ij==0 for j_ij in J_i]):
            return 'no state dependence'
        # an input that depends on state variables in a linear way
        # has a constant derivatives with respect to all state variables
        # (the derivative has no state variables in its free symbols)
        J_ifss=J_i.free_symbols
        svs=set([e for e in self.state_vector])
        inter=J_ifss.intersection(svs)
        if len(inter)==0:
            return 'linear'
        else:
            return 'nonlinear'

    @property
    def no_input_model(self):
        return SmoothReservoirModel(
            self.state_vector,
            self.time_symbol,
            {},# no input fluxes
            self.output_fluxes,
            self.internal_fluxes
        )

    @property
    def function_expressions(self):
        """ Returns the superset of the free symbols of the flux expressions.
        """
        flux_list=self.all_fluxes()
        fun_sets=[ fun_set
                # the sympify in the next line is only necessary for
                # fluxexpressions that are integers (which have no atoms method)
                for fun_set in map(lambda
                    flux:sympify(flux).atoms(Function),flux_list)]
        if len(fun_sets)==0:
            res=set()
        else:
            res=reduce( lambda A,B: A.union(B),fun_sets)
        return res

    @property
    def free_symbols(self):
        """ Returns the superset of the free symbols of the flux expressions including the state variables.
        """
        flux_exprs=self.all_fluxes()
        free_sym_sets=[ sym_set
                # the sympification in the next line is only necessary for
                # fluxexpressions that are numbers
                # It does no harm on expressions
                for sym_set in map(lambda sym:sympify(sym).free_symbols,flux_exprs)]

        if len(free_sym_sets)==0:
            res=set()
        else:
            res=reduce( lambda A,B: A.union(B),free_sym_sets)
        return res

    def subs(self,parameter_dict):
        """ Returns a new instance of class: `SmoothReservoirModel` with all parameters in the parameter_dict replaced
            by their values by calling subs on all the flux expressions.
            Args:
                parameter_dict: A dictionary with the structure {parameter_symbol:parameter_value,....}
        """
        return SmoothReservoirModel(
            self.state_vector,
            self.time_symbol,
            {k:fl.subs(parameter_dict) for k,fl in    self.input_fluxes.items()},
            {k:fl.subs(parameter_dict) for k,fl in   self.output_fluxes.items()},
            {k:fl.subs(parameter_dict) for k,fl in self.internal_fluxes.items()}
        )

    def __str__(self):
        """ This method is called implicitly by print and gives an returns a string that gives an overview over the fluxes
        """
        s = "Object of class "+str(self.__class__)
        indent=2
        s += "\n Input fluxes:\n"
        s += hr.flux_dict_string(self.input_fluxes, indent)

        s += "\n Internal fluxes:\n"
        s += hr.flux_dict_string(self.internal_fluxes, indent)

        s += "\n Output fluxes:\n"
        s += hr.flux_dict_string(self.output_fluxes, indent)

        return s

    def all_fluxes(self):
        # since input and output fluxes are indexed by integers they could
        # overload each other in a common dictionary
        # to avoid this we create a list
        return [v for v in self.input_fluxes.values()] + [v for v in self.output_fluxes.values()] + [v for v in self.internal_fluxes.values()]

    @property
    def jacobian(self):
        state_vec=Matrix(self.state_vector)
        vec=Matrix(self.F)
        return hr.jacobian(vec,state_vec)

    @property
    def is_compartmental(self):
        """ Returns checks that all fluxes are nonnegative
        at the time of implementation this functionality sympy did not support
        relations in predicates yet.
        So while the following works:

        with assuming(Q.positive(x) & Q.positive(y)):
           print(ask(Q.positive(2*x+y)

        it is not possible yet to get a meaningful answer to:

        with assuming(Q.is_true(x>0) & Q.is_true(y>0)):
           print(ask(Q.positive(2*x+y)

        We therefore cannot implement more elaborate assumptions like k_1-(a_12+a_32)>=0
        but still can assume all the state_variables and the time_symbol to be nonnegative
        Therefore we can check the compartmental_property best after all paramater value have been substituted.
        At the moment the function throws an exception if this is not the case.
        """
        #check if all free symbols have been removed
        allowed_symbs= set( [sym for sym in self.state_vector])
        if hasattr(self,"time_symbol"):
            allowed_symbs.add(self.time_symbol)

        if not(allowed_symbs.issuperset(self.free_symbols)):
            raise Exception(
                    Template("Sympy can not check the parameters without assumptions. Try to substitute all variables except the state variables and the time symbol. Use the subs methot of the class {c}").subs(c=self__class__)
                )

        def f(expr):
            res= ask(Q.nonnegative(expr))
            if res is None:
                raise Exception(
                        Template("""Sympy can not (yet) check the parameters even with correct assumptions,\
since relations (<,>) are not implemented yet.
It gave up for the following expression: ${e}."""
                        ).substitute(e=expr)
                    )
            return res

        # making a list of predicated stating that all state variables are nonnegative
        predList=[Q.nonnegative(sym) for sym in self.state_vector]
        if hasattr(self,"time_symbol"):
            predList+=[Q.nonnegative(self.time_symbol)]

        with assuming(*predList):
            # under this assumption eveluate all fluxes
            all_fluxes_nonnegative=all(map(f,self.all_fluxes()))

        return all_fluxes_nonnegative

    # alternative constructor based on the formulation f=u+Bx
    @classmethod
    def from_B_u(cls, state_vector, time_symbol, B, u)->'SmoothReservoirModel':
        """Construct and return a :class:`SmoothReservoirModel` instance from
           :math:`\\dot{x}=B\\,x+u`

        Args:
            state_vector (SymPy dx1-matrix): The model's state vector
                :math:`x`.
                Its entries are SymPy symbols.
            time_symbol (SymPy symbol): The model's time symbol.
            B (SymPy dxd-matrix): The model's compartmental matrix.
            u (SymPy dx1-matrix): The model's external input vector.

        Returns:
            :class:`SmoothReservoirModel`
        """
        #  if not(u):
        #     # fixme mm:
        #     # make sure that ReservoirModels standard constructor can handle an
        #     # empty dict and produce the empty matrix only if necessary later
        #     u=zeros(x.rows,1)
        # fixme mm:
        # we do not seem to have a check that makes sure
        # that the argument B is compartmental
        # maybe the fixme belongs rather to the SmoothModelRun class since
        # we perhaps need parameters

        input_fluxes = hr.in_fluxes_by_index(state_vector, u)
        output_fluxes = hr.out_fluxes_by_index(state_vector, B)
        internal_fluxes = hr.internal_fluxes_by_index(state_vector, B)
        # call the standard constructor
        srm = SmoothReservoirModel(state_vector, time_symbol,
                  input_fluxes, output_fluxes, internal_fluxes)

        return srm

    @property
    def state_variable_set(self):
        return set(self.state_vector)

    @property
    def F(self):
        """SymPy dx1-matrix: The right hand side of the differential equation
        :math:`\\dot{x}=B\\,x+u`."""
        v = (self.external_inputs + self.internal_inputs
                - self.internal_outputs - self.external_outputs)
        #for i in range(len(v)):
        #    v[i] = simplify(v[i])
        return v

    @property
    def external_inputs(self):
        """SymPy dx1-matrix: Return the vector of external inputs."""
        u = zeros(self.nr_pools, 1)
        for k, val in self.input_fluxes.items():
            u[k] = val
        return u

    @property
    def external_outputs(self):
        """SymPy dx1-matrix: Return the vector of external outputs."""
        o = zeros(self.nr_pools, 1)
        for k, val in self.output_fluxes.items():
            o[k] = val
        return o

    @property
    def internal_inputs(self):
        """SymPy dx1-matrix: Return the vector of internal inputs."""
        n = self.nr_pools
        u_int = zeros(n, 1)
        for ln in range(n):
            # find all entries in the fluxes dict that have the target key==ln
            expr = 0
            for k, val in self.internal_fluxes.items():
                if k[1] == ln: #the second part of the tupel is the recipient
                    expr += val
            u_int[ln] = expr
        return u_int

    @property
    def internal_outputs(self):
        """SymPy dx1-matrix: Return the vector of internal outputs."""
        n = self.nr_pools
        o_int = zeros(n, 1)
        for ln in range(n):
            # find all entries in the fluxes dict that have the target key==ln
            expr = 0
            for k, val in self.internal_fluxes.items():
                if k[0] == ln:# the first part of the tupel is the donator
                    expr += val
            o_int[ln] = expr
        return o_int

    @property
    def nr_pools(self):
        """int: Return the number of pools involved in the model."""
        return(len(self.state_variables))

    def port_controlled_Hamiltonian_representation(self):
        """tuple: :math:`J, R, N, x, u` from
        :math:`\\dot{x} = [J(x)-R(x)] \\frac{\\partial}{\\partial x}H+u`.
    	with :math:`H=\\sum_i x_i \\implies \\frac{\\partial}{\\partial x}H =(1,1,...,1)`

        Returns:
            tuple:
            - J (skew symmetric SymPy dxd-matrix) of internal fluxbalances:
              :math:`J_{i,j}=r_{j,i}-r_{i,j}`
            - Q (SymPy dxd-matrix): Diagonal matrix describing the dissipation
              rates (outfluxes).
            - x (SymPy dx1-matrix): The model's state vector.
            - u (SymPy dx1-matrix): The model's external input vector.
        """
        nr_pools = self.nr_pools
        inputs = self.input_fluxes
        outputs = self.output_fluxes
        internal_fluxes = self.internal_fluxes

        C = self.state_vector

        # convert inputs
        u = self.external_inputs

        # calculate decomposition operators
        decomp_fluxes = []
        for pool in range(nr_pools):
            if pool in outputs.keys():
                decomp_flux = outputs[pool]
            else:
                decomp_flux = 0
            decomp_fluxes.append(simplify(decomp_flux))

        Q = diag(*decomp_fluxes)

        # calculate the skewsymmetric matrix J
        J = zeros(nr_pools)

        for (i,j), flux in internal_fluxes.items():
            J[j,i] +=flux
            J[i,j] -=flux

        return (J, Q, C, u)

    def xi_T_N_u_representation(self, factor_out_xi=True):
        """tuple: :math:`\\xi, T, N, x, u` from
        :math:`\\dot{x} = \\xi\\,T\\,N\\,x+u`.

        Args:
            factor_out_xi (bool): If true, xi is extracted from the matrix,
                otherwise :math:`xi=1` will be returned.
                (Defaults to ``True``.)

        Returns:
            tuple:
            - xi (SymPy number): Environmental coefficient.
            - T (SymPy dxd-matrix): Internal fluxes. Main diagonal contains
                ``-1`` entries.
            - N (SymPy dxd-matrix): Diagonal matrix containing the decomposition
                rates.
            - x (SymPy dx1-matrix): The model's state vector.
            - u (SymPy dx1-matrix): The model's external input vector.
        """
        nr_pools = self.nr_pools
        inputs = self.input_fluxes
        outputs = self.output_fluxes
        internal_fluxes = self.internal_fluxes

        C = self.state_vector

        # convert inputs
        u = self.external_inputs

        R = hr.release_operator_1(
            outputs,
            internal_fluxes,
            C
        )
        # calculate transition operator
        T = hr.transfer_operator_3(
            internal_fluxes,
            R,
            C
        )
        # try to extract xi from N and T
        if factor_out_xi:
            xi = hr.factor_out_from_matrix(R)
            # Note mm 02/17/2021
            # since T has -1 on the main diagonal
            # the gcd will be always one so the
            # factor_out_from_matrix(T) is not
            # necessarry.
        else:
            xi = 1

        N = R/xi
        return (xi, T, N, C, u)

    @property
    def compartmental_matrix(self):
        """SymPy Matrix: :math:`B` from
        :math:`\\dot{x} = B\\,x+u`.

        Returns:
            SymPy dxd-matrix: :math:`B = \\xi\\,T\\,N`
        """
        # we could also use the more expensive
        # xi, T, N, C, u = self.xi_T_N_u_representation(factor_out_xi=False))
        return hr.compartmental_matrix_1(
            self.output_fluxes,
            self.internal_fluxes,
            self.state_vector
        )

    def age_moment_system(self, max_order):
        """Return the age moment system of the model.

        Args:
            max_order (int): The maximum order up to which the age moment
                system is created (1 for the mean).

        Returns:
            tuple:
            - extended_state (SymPy d*(max_order+1)x1-matrix): The extended
                state vector of the age moment system.
            - extended_rhs (SymPy d*(max_order+1)x1-matrix): The extended right
                hand side of the age moment ODE.
        """
        u = self.external_inputs
        #X = Matrix(self.state_variables)
        X = self.state_vector
        B = self.compartmental_matrix

        n = self.nr_pools
        extended_state = list(X)
        former_additional_states = [1]*n
        extended_rhs = list(self.F)
        for k in range(1, max_order+1):
            additional_states = [Symbol(str(x)+'_moment_'+str(k)) for x in X]
            g = [k*former_additional_states[i]
                    +(sum([(additional_states[j]-additional_states[i])
                        *B[i,j]*X[j] for j in range(n)])
                      -additional_states[i]*u[i])/X[i] for i in range(n)]

            former_additional_states = additional_states
            extended_state.append(additional_states)
            extended_rhs.append(g)

        extended_state = Matrix(flatten(extended_state))
        extended_rhs = Matrix(flatten(extended_rhs))

        return (extended_state, extended_rhs)

    def plot_pools_and_fluxes(self, ax, mutation_scale = 50, fontsize = 24, thumbnail = False, legend=True, color_fluxes=True, black_and_white=False):
        ax.set_axis_off()
        arrowstyle = "simple"
        visible_pool_names = True

        if color_fluxes:
            pipe_colors = {
                 'linear': 'blue',
                 'nonlinear': 'green',
                 'no state dependence': 'red',
                 'undetermined': 'grey'
            }
        else:
            pipe_colors = {
                 'linear': 'blue',
                 'nonlinear': 'blue',
                 'no state dependence': 'blue',
                 'undetermined': 'blue'
            }

        if thumbnail:
            arrowstyle = "-"
            visible_pool_names = False

        # the former implementation of *_flux_type did not work with Piecewise expressions
        # and had other issues (see the test). To avoid errors stopping the plotting
        # we create a forgiving version of the functions
        def save_maker(func):
            def resilient_func(*args,**kwargs):
                try:
                    res=func(*args,**kwargs)
                except Exception as e:
                    res='undetermined'
                    print(e)
                return res
            return resilient_func

        _res_input_flux_type = save_maker(self._input_flux_type)
        _res_output_flux_type = save_maker(self._output_flux_type)
        _res_internal_flux_type = save_maker(self._internal_flux_type)

        csp = CSPlotter(
            self.state_vector,
            {
                k: _res_input_flux_type(k) for k in self.input_fluxes
                if self.input_fluxes[k] != 0
            },
            {
                k: _res_output_flux_type(k) for k in self.output_fluxes
                if self.output_fluxes[k] != 0
            },
            {
                k: _res_internal_flux_type(*k) for k in self.internal_fluxes
                if self.internal_fluxes[k] != 0
            },
            pipe_colors,
            visible_pool_names = visible_pool_names,
            arrowstyle = arrowstyle,
            fontsize = fontsize
        )
        csp.plot_pools_and_fluxes(ax, black_and_white)

        if legend:
            csp.legend(ax)

    def figure(self, figure_size = (7,7), logo = False, thumbnail = False):
        """Return a figure representing the reservoir model.

        Args:
            figure_size (2-tuple, optional): Width and height of the figure.
                Defaults to (7,7).
            logo (bool, optional): If True, figure_size set to (3,3), no legend,
                smaller font size. Defaults to False.
            thumbnail (bool, optional): If True, produce a very small version,
                no legend. Defaults to False.

        Returns:
            Matplotlib figure: Figure representing the reservoir model.
        """
        fontsize = 24
        mutation_scale = 50
        #mutation_scale=20
        arrowstyle = "simple"
        fontsize = 24
        legend = True

        if thumbnail:
            mutation_scale = 10
            legend = False
            arrowstyle = "-"
            figure_size = (0.7,0.7)

        if logo:
            mutation_scale = 15
            legend = False
            fontsize = 16
            figure_size = (3,3)

        fig = plt.figure(figsize=figure_size, dpi=300)

        if legend:
            #ax = fig.add_axes([0,0,1,0.9])
            ax = fig.add_axes([0,0,0.8,0.8])
        else:
            #ax = fig.add_axes([0,0,1,1])
            ax = fig.add_subplot(1,1,1)

        self.plot_pools_and_fluxes(ax,  mutation_scale = mutation_scale, fontsize = fontsize, thumbnail = thumbnail, legend = legend)

        return fig

def nxgraphs(self):
    return hr.nxgraphs(
        self.state_vector,
        self.inFluxes,
        self.internalFluxes,
        self.outFluxes,
    )

    ##### 14C methods #####

#    def to_14C_only(self, decay_symbol_name, Fa_expr_name):
#        """Construct and return a :class:`SmoothReservoirModel` instance that
#           models the 14C component of the original model.
#
#        Args:
#            decay_symbol_name (str): The name of the 14C decay rate symbol.
#            Fa_expr_name(str): The name of the symbol to be used for the
#                atmospheric C14 fraction function.
#        Returns:
#            :class:`SmoothReservoirModel`
#        """
##        state_vector_14C = Matrix(
##            self.nr_pools,
##            1,
##            [Symbol(sv.name+'_14C') for sv in self.state_vector]
##        )
#        state_vector_14C = self.state_vector
#
#        decay_symbol = Symbol(decay_symbol_name)
#        B_14C = copy(self.compartmental_matrix) - decay_symbol*eye(self.nr_pools)
#        u = self.external_inputs
#        Fa_expr = Function(Fa_expr_name)(self.time_symbol)
#        u_14C = Matrix(self.nr_pools, 1, [expr*Fa_expr for expr in u])
#
#        srm_14C = SmoothReservoirModel.from_B_u(
#            state_vector_14C,
#            self.time_symbol,
#            B_14C,
#            u_14C
#        )
#
#        return srm_14C
#
#    def to_14C_explicit(self, decay_symbol_name, Fa_expr_name):
#        """Construct and return a :class:`SmoothReservoirModel` instance that
#           models the 14C component additional to the original model.
#
#        Args:
#            decay_symbol_name (str): The name of the 14C decay rate symbol.
#            Fa_expr_name(str): The name of the symbol to be used for the
#                atmospheric C14 fraction function.
#        Returns:
#            :class:`SmoothReservoirModel`
#        """
#        state_vector = self.state_vector
#        B, u = self.compartmental_matrix, self.external_inputs
#        srm_14C = self.to_14C_only(decay_symbol_name, Fa_expr_name)
#        state_vector_14C = srm_14C.state_vector
#        B_C14 = srm_14C.compartmental_matrix
#        u_14C = srm_14C.external_inputs
#
#        nr_pools = self.nr_pools
#
#        state_vector_total = Matrix(nr_pools*2, 1, [1]*(nr_pools*2))
#        state_vector_total[:nr_pools,0] = state_vector
#        state_vector_total[nr_pools:,0] = state_vector_14C
#
#        B_total = eye(nr_pools*2)
#        B_total[:nr_pools, :nr_pools] = B
#        B_total[nr_pools:, nr_pools:] = B_C14
#
#        u_total = Matrix(nr_pools*2, 1, [1]*(nr_pools*2))
#        u_total[:nr_pools,0] = u
#        u_total[nr_pools:,0] = u_14C
#
#        srm_total = SmoothReservoirModel.from_B_u(
#            state_vector_total,
#            self.time_symbol,
#            B_total,
#            u_total)
#
#        return srm_total

    def steady_states(self, par_set = None):
        if par_set is None:
            #compute steady state formulas
            par_set = {}
        # try to calculate the steady states for ten seconds
        # after ten seconds stop it
        q = multiprocessing.Queue()
        def calc_steady_states(q):
            ss = solve(self.F.subs(par_set), self.state_vector, dict=True)
            q.put(ss)

        p = multiprocessing.Process(target=calc_steady_states, args=(q,))
        p.start()
        p.join(10)
        if p.is_alive():
            p.terminate()
            p.join()
            steady_states = []
        else:
            steady_states = q.get()

        formal_steady_states = []
        for ss in steady_states:
            result = []
            ss_dict = {}
            for sv_symbol in self.state_vector:
                if sv_symbol in ss.keys():
                    ss[sv_symbol] = simplify(ss[sv_symbol])
                else:
                    ss[sv_symbol] = sv_symbol

                ss_expr = ss[sv_symbol]
                if self.time_symbol in ss_expr.free_symbols:
                    # take limit of time to infinity if steady state still depends on time
                    ss_expr = limit(ss_expr, self.time_symbol, oo)
                ss_dict[sv_symbol.name] = ss_expr

            formal_steady_states.append(ss_dict)

        return formal_steady_states

    ##### functions for internal use only #####

max_stars_repo_name: goujou/CompartmentalSystems

max_issues_repo_name: goujou/CompartmentalSystems

max_forks_repo_name: goujou/CompartmentalSystems

reason: import numpy,from sympy,import networkx


text:
from methylprep.processing import pipeline
import pandas as pd
from pathlib import Path
import numpy as np

class TestBatchSize():

    """ TOO SLOW for CI - and redundant with test_pipeline
    @staticmethod
    def test_no_batch_size():
        test_data_dir = 'docs/example_data/GSE69852'
        df = pipeline.run_pipeline(test_data_dir, export=True)
        if not np.isclose(df[0].unmethylated.data_frame.iloc[0]['noob'], 4119.633578946326):
            raise AssertionError()
        if not Path('docs/example_data/GSE69852/9247377093/9247377093_R02C01_processed.csv').is_file():
            raise AssertionError()

    @staticmethod
    def test_no_batch_size_betas():
        test_data_dir = 'docs/example_data/GSE69852'
        betas = pipeline.run_pipeline(test_data_dir, betas=True)
        if not np.isclose(betas.iloc[0]['9247377093_R02C01'], 0.23623395577166542):
            raise AssertionError()
        if not Path('beta_values.pkl').is_file():
            raise AssertionError()

    @staticmethod
    def test_no_batch_size_m_values():
        test_data_dir = 'docs/example_data/GSE69852'
        m_values = pipeline.run_pipeline(test_data_dir, m_value=True)
        if not np.isclose(m_values.iloc[0]['9247377093_R02C01'], -1.6575579144514734):
            raise AssertionError()
        if not Path('m_values.pkl').is_file():
            raise AssertionError()
    """

    @staticmethod
    def test_with_batch_size():
        test_data_dir = 'docs/example_data/GSE69852'
        df = pipeline.run_pipeline(test_data_dir, export=True, batch_size=1, sample_name='AdultLiver1')
        #if not np.isclose(df[0].unmethylated.data_frame.iloc[0]['noob'], 4119.633578946326):
        if not np.isclose(df[0]._SampleDataContainer__data_frame.iloc[0]['beta_value'], 0.236):
            raise AssertionError()
        if not np.isclose(df[0]._SampleDataContainer__data_frame.iloc[2]['m_value'], 4.146):
            raise AssertionError()
        if not Path('docs/example_data/GSE69852/9247377093/9247377093_R02C01_processed.csv').is_file():
            raise AssertionError()
        print(f"TEST OUTPUT FILES: {list(Path(test_data_dir).rglob('*'))}")
        for file in Path(test_data_dir).rglob('*.pkl'):
            file.unlink()

    @staticmethod
    def test_batch_size_betas():
        test_data_dir = 'docs/example_data/GSE69852'
        betas = pipeline.run_pipeline(test_data_dir, betas=True, batch_size=1)
        if not np.isclose(betas.iloc[0]['9247377093_R02C01'], 0.23624517): #0.23623395577166542):
            raise AssertionError(f"{betas.iloc[0]['9247377093_R02C01']} != 0.23623395577166542")
        if not Path(test_data_dir, 'beta_values.pkl').is_file():
            raise AssertionError()
        print(f"TEST OUTPUT FILES: {list(Path(test_data_dir).rglob('*'))}")
        for file in Path(test_data_dir).rglob('*.pkl'):
            file.unlink()

    """
    @staticmethod
    def test_batch_size_m_values():
        test_data_dir = 'docs/example_data/GSE69852'
        m_values = pipeline.run_pipeline(test_data_dir, m_value=True, batch_size=1)
        if not np.isclose(m_values.iloc[0]['9247377093_R02C01'], -1.6575579144514734):
            raise AssertionError()
        if not (Path('m_values_1.pkl').is_file() and Path('m_values_2.pkl').is_file()):
            raise AssertionError()
        if not Path('m_values.pkl').is_file():
            raise AssertionError()
    """

ext: py

max_stars_repo_name: holonomicjl/methylprep

max_issues_repo_name: holonomicjl/methylprep

max_forks_repo_name: holonomicjl/methylprep

reason: import numpy


text:
import cv2
from sklearn import svm
import os
import numpy as np
from sklearn.externals import joblib
from skimage.feature import hog
import argparse

parser = argparse.ArgumentParser(description='Parse Training Directory')
parser.add_argument('--pos', help='Path to directory containing Positive Images')
parser.add_argument('--neg', help='Path to directory containing Negative images')

args = parser.parse_args()

pos_img_dir = args.pos
neg_img_dir = args.neg

clf = joblib.load('person_final.pkl')

total_pos_samples = 0
total_neg_samples = 0

def crop_centre(img):
    h, w, d = img.shape
    l = (w - 64)/2
    t = (h - 128)/2
    #print (h, w, l, t)
    crop = img[t:t+128, l:l+64]
    return crop

def read_filenames():

    f_pos = []
    f_neg = []

    for (dirpath, dirnames, filenames) in os.walk(pos_img_dir):
        f_pos.extend(filenames)
        break

    for (dirpath, dirnames, filenames) in os.walk(neg_img_dir):
        f_neg.extend(filenames)
        break

    print "Positive Image Samples: " + str(len(f_pos))
    print "Negative Image Samples: " + str(len(f_neg))

    return f_pos, f_neg

def read_images(f_pos, f_neg):

    print ("Reading Images")

    array_pos_features = []
    array_neg_features = []
    global total_pos_samples
    global total_neg_samples
    for imgfile in f_pos:
        img = cv2.imread(os.path.join(pos_img_dir, imgfile))
        cropped = crop_centre(img)
        gray = cv2.cvtColor(cropped, cv2.COLOR_BGR2GRAY)
        features = hog(gray, orientations=9, pixels_per_cell=(8, 8), cells_per_block=(2, 2), block_norm="L2", feature_vector=True)
        array_pos_features.append(features.tolist())

        total_pos_samples += 1

    for imgfile in f_neg:
        img = cv2.imread(os.path.join(neg_img_dir, imgfile))
        cropped = crop_centre(img)
        gray = cv2.cvtColor(cropped, cv2.COLOR_BGR2GRAY)
        features = hog(gray, orientations=9, pixels_per_cell=(8, 8), cells_per_block=(2, 2), block_norm="L2", feature_vector=True)
        array_neg_features.append(features.tolist())
        total_neg_samples += 1

    return array_pos_features, array_neg_features

pos_img_files, neg_img_files = read_filenames()

pos_features, neg_features = read_images(pos_img_files, neg_img_files)

pos_result = clf.predict(pos_features)
neg_result = clf.predict(neg_features)

true_positives = cv2.countNonZero(pos_result)
false_negatives = pos_result.shape[0] - true_positives

false_positives = cv2.countNonZero(neg_result)
true_negatives = neg_result.shape[0] - false_positives

print ("True Positives: " + str(true_positives), "False Positives: " + str(false_positives))
print ("True Negatives: " + str(true_negatives), "False Negatives: " + str(false_negatives))

precision = float(true_positives) / (true_positives + false_positives)
recall = float(true_positives) / (true_positives + false_negatives)

f1 = 2*precision*recall / (precision + recall)

print ("Precision: " + str(precision), "Recall: " + str(recall))
print ("F1 Score: " + str(f1))

max_stars_repo_name: enesdoruk/Autonom-Vehicles-Projects

max_issues_repo_name: enesdoruk/Autonom-Vehicles-Projects

max_forks_repo_name: enesdoruk/Autonom-Vehicles-Projects


text:
'''
// NACA 2412 Airfoil
// Created using Gmsh (http://gmsh.info/)
// Mesh Distributed by: Notre Dame CICS (MIT Liscense)
// - Associated publication:
// doi: https://doi.org/10.1016/j.jcp.2019.01.021
// github: https://github.com/cics-nd/rans-uncertainty
'''

import pygmsh #See: https://pypi.org/project/pygmsh/
import numpy as np
import meshio

def createMesh():
    # Airfoil coordinates
    airfoil_coordinates = np.genfromtxt('airfoil-profile.data', delimiter=',', skip_header=1, skip_footer=0)
    airfoil_coordinates = np.c_[airfoil_coordinates, np.zeros(airfoil_coordinates.shape[0])]
    print(airfoil_coordinates)
    # Scale airfoil to input coord
    coord = 1.0
    airfoil_coordinates *= coord

    # Instantiate geometry object
    geom = pygmsh.built_in.Geometry()

    # Create polygon for airfoil
    char_length = 5.0e-2 # Resolution
    airfoil = geom.add_polygon(airfoil_coordinates, char_length, make_surface=False)

    # Create surface for numerical domain with an airfoil-shaped hole
    left_dist = 1.0
    right_dist = 3.0
    top_dist = 1.0
    bottom_dist = 1.0
    xmin = airfoil_coordinates[:, 0].min() - left_dist * coord
    xmax = airfoil_coordinates[:, 0].max() + right_dist * coord
    ymin = airfoil_coordinates[:, 1].min() - bottom_dist * coord
    ymax = airfoil_coordinates[:, 1].max() + top_dist * coord
    domainCoordinates = np.array(
        [[xmin, ymin, 0.0], [xmax, ymin, 0.0], [xmax, ymax, 0.0], [xmin, ymax, 0.0]]
    )
    polygon = geom.add_polygon(domainCoordinates, char_length, holes=[airfoil])

    geom.add_raw_code("Recombine Surface {{{}}};".format(polygon.surface.id))

    axis = [0, 0, 1]
    geom.extrude(
        polygon,
        translation_axis=axis,
        num_layers=3.0/char_length,
        recombine=True
    )

    ref = 10.525891646546
    points, cells, _, _, _ = pygmsh.generate_mesh(geom, remove_faces=True, gmsh_path='../gmsh', geo_filename='LES-mesh.geo')
    return points, cells

if __name__ == "__main__":

    createMesh()
    # meshio.write_points_cells("airfoil.vtu", *test())

max_issues_repo_name: zabaras/rans-uncertainty

max_forks_repo_name: zabaras/rans-uncertainty


text:
#!/usr/bin/env python
import rospy
import numpy as np
import atexit
import tf
from os.path import expanduser
from time import gmtime, strftime
from numpy import linalg as LA
from tf.transformations import euler_from_quaternion
from nav_msgs.msg import Odometry

home = expanduser('~')
file = open(strftime(home+'/rcws/logs/wp-%Y-%m-%d-%H-%M-%S',gmtime())+'.csv', 'w')

def save_waypoint(data):
    quaternion = np.array([data.pose.pose.orientation.x,
                           data.pose.pose.orientation.y,
                           data.pose.pose.orientation.z,
                           data.pose.pose.orientation.w])

    euler = tf.transformations.euler_from_quaternion(quaternion)
    speed = LA.norm(np.array([data.twist.twist.linear.x,
                              data.twist.twist.linear.y,
                              data.twist.twist.linear.z]),2)
    if data.twist.twist.linear.x>0.:
        print data.twist.twist.linear.x

    file.write('%f, %f, %f, %f\n' % (data.pose.pose.position.x,
                                     data.pose.pose.position.y,
                                     euler[2],
                                     speed))

def shutdown():
    file.close()
    print('Goodbye')

def listener():
    rospy.init_node('waypoints_logger', anonymous=True)
    rospy.Subscriber('pf/pose/odom', Odometry, save_waypoint)
    rospy.spin()

if __name__ == '__main__':
    atexit.register(shutdown)
    print('Saving waypoints...')
    try:
        listener()
    except rospy.ROSInterruptException:
        pass

reason: import numpy,from numpy


text:
import numpy as np
from tensorflow import keras

def is_numpy(obj):
    """
    Check of the type is instance of numpy array
    :param obj: object to check
    :return: True if the object is numpy-type array.
    """
    return isinstance(obj, (np.ndarray, np.generic))

def ensure_numpy_type(obj):
    """
    Raise exception if it's not a numpy
    :param obj: object to check
    :return: numpy object
    """
    if is_numpy(obj):
        return obj
    else:
        raise AttributeError('Not a numpy type.')

def ensure_tf_type(obj, fake_input_layer=None, name=None):
    """
    Convert to Keras Constant if needed
    :param obj: numpy / tf type
    :param fake_input_layer: fake input layer to add constant
    :return: tf type
    """
    if is_numpy(obj):
        if obj.dtype == np.int64:
            obj = np.int32(obj)

        def target_layer(_, inp=obj, dtype=obj.dtype.name):
            import numpy as np
            import tensorflow as tf
            if not isinstance(inp, (np.ndarray, np.generic)):
                inp = np.array(inp, dtype=dtype)
            return tf.constant(inp, dtype=inp.dtype)

        lambda_layer = keras.layers.Lambda(target_layer, name=name)
        return lambda_layer(fake_input_layer)
    else:
        return obj

def check_torch_keras_error(model, k_model, input_np, epsilon=1e-5, change_ordering=False):
    """
    Check difference between Torch and Keras models
    :param model: torch model
    :param k_model: keras model
    :param input_np: input data as numpy array or list of numpy array
    :param epsilon: allowed difference
    :param change_ordering: change ordering for keras input
    :return: actual difference
    """
    from torch.autograd import Variable
    import torch

    if isinstance(input_np, np.ndarray):
        input_np = [input_np]

    input_var = [Variable(torch.FloatTensor(i)) for i in input_np]
    pytorch_output = model(*input_var)
    if not isinstance(pytorch_output, tuple):
        pytorch_output = [pytorch_output.data.numpy()]
    else:
        pytorch_output = [p.data.numpy() for p in pytorch_output]

    if change_ordering:
        # change image data format
        _input_np = []
        for i in input_np:
            axes = list(range(len(i.shape)))
            axes = axes[0:1] + axes[2:] + axes[1:2]
            _input_np.append(np.transpose(i, axes))
        input_np = _input_np

        # run keras model
        keras_output = k_model.predict(input_np)
        if not isinstance(keras_output, list):
            keras_output = [keras_output]

        # change image data format
        _koutput = []
        for k in keras_output:
            axes = list(range(len(k.shape)))
            axes = axes[0:1] + axes[-1:] + axes[1:-1]
            _koutput.append(np.transpose(k, axes))
        keras_output = _koutput
    else:
        keras_output = k_model.predict(input_np)
        if not isinstance(keras_output, list):
            keras_output = [keras_output]

    max_error = 0
    for p, k in zip(pytorch_output, keras_output):
        error = np.max(np.abs(p - k))
        if error > max_error:
            max_error = error

    assert max_error < epsilon
    return max_error


text:
import re, os, ntpath
import numpy as np
from . import utils

channelmap = {
    'Xrotation': 'x',
    'Yrotation': 'y',
    'Zrotation': 'z'
}

channelmap_inv = {
    'x': 'Xrotation',
    'y': 'Yrotation',
    'z': 'Zrotation',
}

ordermap = {
    'x': 0,
    'y': 1,
    'z': 2,
}

class Anim(object):
    """
    A very basic animation object
    """
    def __init__(self, quats, pos, offsets, parents, bones):
        """
        :param quats: local quaternions tensor
        :param pos: local positions tensor
        :param offsets: local joint offsets
        :param parents: bone hierarchy
        :param bones: bone names
        """
        self.quats = quats
        self.pos = pos
        self.offsets = offsets
        self.parents = parents
        self.bones = bones

def read_bvh(filename, start=None, end=None, order=None):
    """
    Reads a BVH file and extracts animation information.

    :param filename: BVh filename
    :param start: start frame
    :param end: end frame
    :param order: order of euler rotations
    :return: A simple Anim object conatining the extracted information.
    """

    f = open(filename, "r")

    i = 0
    active = -1
    end_site = False

    names = []
    orients = np.array([]).reshape((0, 4))
    offsets = np.array([]).reshape((0, 3))
    parents = np.array([], dtype=int)

    # Parse the  file, line by line
    for line in f:

        if "HIERARCHY" in line: continue
        if "MOTION" in line: continue

        rmatch = re.match(r"ROOT (\w+)", line)
        if rmatch:
            names.append(rmatch.group(1))
            offsets = np.append(offsets, np.array([[0, 0, 0]]), axis=0)
            orients = np.append(orients, np.array([[1, 0, 0, 0]]), axis=0)
            parents = np.append(parents, active)
            active = (len(parents) - 1)
            continue

        if "{" in line: continue

        if "}" in line:
            if end_site:
                end_site = False
            else:
                active = parents[active]
            continue

        offmatch = re.match(r"\s*OFFSET\s+([\-\d\.e]+)\s+([\-\d\.e]+)\s+([\-\d\.e]+)", line)
        if offmatch:
            if not end_site:
                offsets[active] = np.array([list(map(float, offmatch.groups()))])
            continue

        chanmatch = re.match(r"\s*CHANNELS\s+(\d+)", line)
        if chanmatch:
            channels = int(chanmatch.group(1))
            if order is None:
                channelis = 0 if channels == 3 else 3
                channelie = 3 if channels == 3 else 6
                parts = line.split()[2 + channelis:2 + channelie]
                if any([p not in channelmap for p in parts]):
                    continue
                order = "".join([channelmap[p] for p in parts])
            continue

        jmatch = re.match("\s*JOINT\s+(\w+)", line)
        if jmatch:
            names.append(jmatch.group(1))
            offsets = np.append(offsets, np.array([[0, 0, 0]]), axis=0)
            orients = np.append(orients, np.array([[1, 0, 0, 0]]), axis=0)
            parents = np.append(parents, active)
            active = (len(parents) - 1)
            continue

        if "End Site" in line:
            end_site = True
            continue

        fmatch = re.match("\s*Frames:\s+(\d+)", line)
        if fmatch:
            if start and end:
                fnum = (end - start) - 1
            else:
                fnum = int(fmatch.group(1))
            positions = offsets[np.newaxis].repeat(fnum, axis=0)
            rotations = np.zeros((fnum, len(orients), 3))
            continue

        fmatch = re.match("\s*Frame Time:\s+([\d\.]+)", line)
        if fmatch:
            frametime = float(fmatch.group(1))
            continue

        if (start and end) and (i < start or i >= end - 1):
            i += 1
            continue

        dmatch = line.strip().split(' ')
        if dmatch:
            data_block = np.array(list(map(float, dmatch)))
            N = len(parents)
            fi = i - start if start else i
            if channels == 3:
                positions[fi, 0:1] = data_block[0:3]
                rotations[fi, :] = data_block[3:].reshape(N, 3)
            elif channels == 6:
                data_block = data_block.reshape(N, 6)
                positions[fi, :] = data_block[:, 0:3]
                rotations[fi, :] = data_block[:, 3:6]
            elif channels == 9:
                positions[fi, 0] = data_block[0:3]
                data_block = data_block[3:].reshape(N - 1, 9)
                rotations[fi, 1:] = data_block[:, 3:6]
                positions[fi, 1:] += data_block[:, 0:3] * data_block[:, 6:9]
            else:
                raise Exception("Too many channels! %i" % channels)

            i += 1

    f.close()

    rotations = utils.euler_to_quat(np.radians(rotations), order=order)
    rotations = utils.remove_quat_discontinuities(rotations)

    return Anim(rotations, positions, offsets, parents, names)

def get_lafan1_set(bvh_path, actors, window=50, offset=20, debug=False):
    """
    Extract the same test set as in the article, given the location of the BVH files.

    :param bvh_path: Path to the dataset BVH files
    :param list: actor prefixes to use in set
    :param window: width  of the sliding windows (in timesteps)
    :param offset: offset between windows (in timesteps)
    :return: tuple:
        X: local positions
        Q: local quaternions
        parents: list of parent indices defining the bone hierarchy
        contacts_l: binary tensor of left-foot contacts of shape (Batchsize, Timesteps, 2)
        contacts_r: binary tensor of right-foot contacts of shape (Batchsize, Timesteps, 2)
    """
    npast = 10
    subjects = []
    seq_names = []
    X = []
    Q = []
    contacts_l = []
    contacts_r = []

    # Extract
    bvh_files = os.listdir(bvh_path)
    if debug:
        bvh_files = bvh_files[:3]

    for file in bvh_files:
        if file.endswith('.bvh'):
            # basename = ntpath.basename(file)[:-4].split('_')
            # seq_name, subject = basename[0], basename[1]
            # seq_name, subject = ntpath.basename(file[:-4]).split('_')
            seq_name = subject = ntpath.basename(file[:-4])

            if subject in actors:
                print('Processing file {}'.format(file))
                seq_path = os.path.join(bvh_path, file)
                anim = read_bvh(seq_path)

                # Sliding windows
                i = 0
                while i+window < anim.pos.shape[0]:
                    q, x = utils.quat_fk(anim.quats[i: i+window], anim.pos[i: i+window], anim.parents)
                    # Extract contacts
                    c_l, c_r = utils.extract_feet_contacts(x, [3, 4], [7, 8], velfactor=0.02)
                    X.append(anim.pos[i: i+window])
                    Q.append(anim.quats[i: i+window])
                    seq_names.append(seq_name)
                    subjects.append(subjects)
                    contacts_l.append(c_l)
                    contacts_r.append(c_r)

                    i += offset

    X = np.asarray(X).astype('float32')
    Q = np.asarray(Q).astype('float32')
    contacts_l = np.asarray(contacts_l).astype('float32')
    contacts_r = np.asarray(contacts_r).astype('float32')

    # Sequences around XZ = 0
    xzs = np.mean(X[:, :, 0, ::2], axis=1, keepdims=True)
    X[:, :, 0, 0] = X[:, :, 0, 0] - xzs[..., 0]
    X[:, :, 0, 2] = X[:, :, 0, 2] - xzs[..., 1]

    # Unify facing on last seed frame
    # X, Q = utils.rotate_at_frame(X, Q, anim.parents, n_past=npast)

    return X.astype('float32'), Q.astype('float32'), anim.parents, contacts_l, contacts_r

def get_train_stats(bvh_folder, train_set):
    """
    Extract the same training set as in the paper in order to compute the normalizing statistics
    :return: Tuple of (local position mean vector, local position standard deviation vector, local joint offsets tensor)
    """
    print('Building the train set...')
    xtrain, qtrain, parents, _, _ = get_lafan1_set(bvh_folder, train_set, window=50, offset=20)

    print('Computing stats...\n')
    # Joint offsets : are constant, so just take the first frame:
    offsets = xtrain[0:1, 0:1, 1:, :]  # Shape : (1, 1, J, 3)

    # Global representation:
    q_glbl, x_glbl = utils.quat_fk(qtrain, xtrain, parents)

    # Global positions stats:
    x_mean = np.mean(x_glbl.reshape([x_glbl.shape[0], x_glbl.shape[1], -1]).transpose([0, 2, 1]), axis=(0, 2), keepdims=True)
    x_std = np.std(x_glbl.reshape([x_glbl.shape[0], x_glbl.shape[1], -1]).transpose([0, 2, 1]), axis=(0, 2), keepdims=True)

    return x_mean, x_std, offsets

max_stars_repo_name: Garfield-kh/PoseTriplet

max_issues_repo_name: Garfield-kh/PoseTriplet

max_forks_repo_name: Garfield-kh/PoseTriplet


text:
from .._C_netket.operator import LocalOperator as _LocalOperator

import numpy as _np

def sigmax(hilbert, site):
    """
    Builds the sigma_x operator acting on the `site`-th of the Hilbert
    space `hilbert`.

    If `hilbert` is a non-Spin space of local dimension M, it is considered
    as a (M-1)/2 - spin space.

    :param hilbert: The hilbert space
    :param site: the site on which this operator acts
    :return: a nk.operator.LocalOperator
    """
    N = hilbert.local_size
    S = (N - 1) / 2

    D = [_np.sqrt((S + 1) * 2 * a - a * (a + 1)) for a in _np.arange(1, N)]
    mat = _np.diag(D, 1) + _np.diag(D, -1)
    return _LocalOperator(hilbert, mat, [site])

def sigmay(hilbert, site):
    """
    Builds the sigma_y operator acting on the `site`-th of the Hilbert
    space `hilbert`.

    If `hilbert` is a non-Spin space of local dimension M, it is considered
    as a (M-1)/2 - spin space.

    :param hilbert: The hilbert space
    :param site: the site on which this operator acts
    :return: a nk.operator.LocalOperator
    """
    N = hilbert.local_size
    S = (N - 1) / 2

    D = _np.array(
        [1j * _np.sqrt((S + 1) * 2 * a - a * (a + 1)) for a in _np.arange(1, N)]
    )
    mat = _np.diag(D, -1) + _np.diag(-D, 1)
    return _LocalOperator(hilbert, mat, [site])

def sigmaz(hilbert, site):
    """
    Builds the sigma_z operator acting on the `site`-th of the Hilbert
    space `hilbert`.

    If `hilbert` is a non-Spin space of local dimension M, it is considered
    as a (M-1)/2 - spin space.

    :param hilbert: The hilbert space
    :param site: the site on which this operator acts
    :return: a nk.operator.LocalOperator
    """
    N = hilbert.local_size
    S = (N - 1) / 2

    D = _np.array([2 * m for m in _np.arange(S, -(S + 1), -1)])
    mat = _np.diag(D, 0)
    return _LocalOperator(hilbert, mat, [site])

def sigmam(hilbert, site):
    """
    Builds the $sigma_- = sigma_x - im * sigma_y$ operator acting on the
    `site`-th of the Hilbert space `hilbert`.

    If `hilbert` is a non-Spin space of local dimension M, it is considered
    as a (M-1)/2 - spin space.

    :param hilbert: The hilbert space
    :param site: the site on which this operator acts
    :return: a nk.operator.LocalOperator
    """
    N = hilbert.local_size
    S = (N - 1) / 2

    S2 = (S + 1) * S
    D = _np.array([_np.sqrt(S2 - m * (m - 1)) for m in _np.arange(S, -S, -1)])
    mat = _np.diag(D, -1)
    return _LocalOperator(hilbert, mat, [site])

def sigmap(hilbert, site):
    """
    Builds the $sigma_- = sigma_x + im * sigma_y$ operator acting on the
    `site`-th of the Hilbert space `hilbert`.

    If `hilbert` is a non-Spin space of local dimension M, it is considered
    as a (M-1)/2 - spin space.

    :param hilbert: The hilbert space
    :param site: the site on which this operator acts
    :return: a nk.operator.LocalOperator
    """
    N = hilbert.local_size
    S = (N - 1) / 2

    S2 = (S + 1) * S
    D = _np.array([_np.sqrt(S2 - m * (m + 1)) for m in _np.arange(S - 1, -(S + 1), -1)])
    mat = _np.diag(D, 1)
    return _LocalOperator(hilbert, mat, [site])

max_stars_repo_name: vigsterkr/netket

max_issues_repo_name: vigsterkr/netket

max_forks_repo_name: vigsterkr/netket


text:
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 15 10:09:03 2021

@author: Rin Sun
"""
import torch
import torch.nn as nn
import torch.nn.functional as F

import time # record the time record
import matplotlib.pyplot as plt #draw the diagram

import numpy as np #linear algebra
import datetime
import random
from math import *
from dataloader import *
from config import *
import config

class Parser(nn.Module):
    def __init__(self,op_dim = 64, args_dim = 32,state_dim = 128, vocab = 3000):
        super(Parser, self).__init__()

        # create a list of vectors representing operators and words
        self.wordvecs = nn.Embedding(vocab, 64)
        self.op_vecs = nn.Embedding(vocab, op_dim)
        self.arg_vecs = nn.Embedding(vocab, args_dim)

        # Input of this model will be a set of embeddings of the operator

        # Recreate operators using pytorch
        self.state_dim = state_dim
        self.args_dim = args_dim
        self.op_dim = op_dim

        # Synthesis methods required
        self.log_p = 0
        self.count = 0
        self.program = ""

        # to create the probability analyzer unit
        self.pfc1 = nn.Linear(op_dim+state_dim,300)
        self.pfc2 = nn.Linear(300,200)
        self.pfca = nn.Linear(200,200)
        self.pfc3 = nn.Linear(200,1)
        # to create the repeater unit
        self.rfc1 = nn.Linear(args_dim+state_dim,200)
        self.rfc2 = nn.Linear(200,399)
        self.rfc3 = nn.Linear(399,state_dim)
        # Gated recurrential unit
        self.GRU = torch.nn.GRU(64,state_dim,1)
        self.namo = "Melian"

    def pdf_ops(self,opp,state):
        # shape of input variables are [num, embed_dim]
        # output will be the pdf tensor and the argmax position
        states = torch.broadcast_to(state,[len(opp),self.state_dim])

        r = torch.cat([opp,states],1)

        pdf_tensor = F.tanh(self.pfc1(r))
        pdf_tensor = F.softplus(self.pfc2(pdf_tensor))
        pdf_tensor = F.softplus(self.pfc3(pdf_tensor))
        #create pdf tensor
        pdf = pdf_tensor/torch.sum(pdf_tensor,0)

        #index of the maximum operator is
        index = np.argmax(pdf.detach().numpy())
        return pdf, index

    def repeat_vector(self,semantics,arg):
        # Create next semantics vector for the argument
        r = torch.cat([semantics,arg],-1)
        r = F.tanh(self.rfc1(r))
        r = F.softplus(self.rfc2(r))
        r = F.tanh(self.rfc3(r))

        return r

    def convert(self,x):
        # Start the parsing process
        self.log_p = 0
        self.program = ""
        def parse(s,arg,ops,ops_dict):
            # analogue goal-based policy in multi-goal RL
            # estimate the action to take when given the state (s) and the goal (arg)
            parse_state = self.repeat_vector(s,arg)
            # 1. create a state given the current state and the goal given
            pdf,index = self.pdf_ops(ops,parse_state)
            self.log_p  = self.log_p - torch.log(pdf[index])
            operator = ops_dict[index+1]
            self.program += str(operator)
            # 2. pass the semantics down to next level if required

            if operator in arged_ops:
                self.program +=  "("

                args = arg_dict[operator]
                args_paramed = self.arg_vecs(torch.tensor(args))

                for i in range(len(args)):

                    parse(parse_state,torch.reshape(args_paramed[i],[1,self.args_dim]),ops,ops_dict)
                    if i != (len(args)-1):
                        self.program += ","
                self.program += ")"

        root = self.arg_vecs(torch.tensor([0]))
        ops = self.op_vecs(torch.tensor(config.ops))

        # TODO: dynamic library useage during the parsing
        parse(x,root,ops,operators)
        return self.program,self.log_p

    def program_prob(self,x,ops_sequence):
        # Start the parsing process
        self.log_p = 0
        self.count = 0
        self.program = ""
        def parse(s,arg,ops,ops_dict,ops_sequence):
            # analogue goal-based policy in multi-goal RL
            # estimate the action to take when given the state (s) and the goal (arg)
            parse_state = self.repeat_vector(s,arg)
            # 1. create a state given the current state and the goal given
            pdf,index = self.pdf_ops(ops,parse_state)
            index = ops_dict.index(ops_sequence[self.count])
            self.log_p  = self.log_p - torch.log(pdf[index-1])
            operator = ops_dict[index]
            self.program += str(operator)
            # 2. pass the semantics down to next level if required
            self.count += 1

            if operator in arged_ops:
                self.program +=  "("

                args = arg_dict[operator]
                args_paramed = self.arg_vecs(torch.tensor(args))

                for i in range(len(args)):

                    parse(parse_state,torch.reshape(args_paramed[i],[1,self.args_dim]),ops,ops_dict,ops_sequence)
                    if i != (len(args)-1):
                        self.program += ","
                self.program += ")"

        root = self.arg_vecs(torch.tensor([0]))
        ops = self.op_vecs(torch.tensor(config.ops))

        # TODO: dynamic library useage during the parsing
        parse(x,root,ops,operators,ops_sequence)
        return self.program,self.log_p

    def save_model(self):
        dirc = "D:/AnQira/" + str(self.namo) + "/"
        torch.save(self, dirc+self.namo + '.pth')

    def load_model(self):
        dirc = "D:/AnQira/" + str(self.namo) + "/"
        try:
            self = torch.load(dirc+self.namo + '.pth')
        except:
            print("Failed")

    def run(self,x,seq = None):
        sequence = self.wordvecs(torch.tensor(x))

        results,hidden  = self.GRU(sequence)

        hidden = torch.reshape(hidden,[hidden.shape[1],128])
        semantics = torch.reshape(hidden[hidden.shape[0]-1],[1,128])

        if seq != None:
            program,loss  = self.program_prob(semantics,seq)
        else:
            program,loss = self.convert(semantics)
        return program, loss

ps = Parser() #Parser that coverts the sequence

Adam = torch.optim.Adam(ps.parameters(),0.0003) # An adam optimizer is initiated

for epoch in range(EPOCH):
    Loss = 0
    for i in range(len(tasks)):
        toks = data.convert_to_data(tasks[i][0])
        program, loss = ps.run(toks,Decompose(tasks[i][1]))

        Loss = Loss + loss
        #print(tasks[i][0],tasks[i][1])
    ps.zero_grad()
    Loss.backward()
    Adam.step()

    print("EPOCH:",epoch,"Loss: ",Loss.detach().numpy())

ps.save_model()
us = Parser()
us.load_model()
num = 0
for i in range(len(tasks)):
    toks = data.convert_to_data(tasks[i][0])
    program, loss = ps.run(toks,Decompose(tasks[i][1]))
    l = loss.detach().numpy()
    print("prob:",exp(-l))
    if (program == tasks[i][1]):
        num = num + 1
    else:
        print(tasks[i][0])
        print(program,tasks[i][1])

def convert_to_program(x,tar):
    toks = data.convert_to_data(x)
    program, loss = tar.run(toks,Decompose(tasks[i][1]))
    return program,exp(-loss)

max_stars_repo_name: Rin-The-QT-Bunny/tree-structrued-parser

max_issues_repo_name: Rin-The-QT-Bunny/tree-structrued-parser

max_forks_repo_name: Rin-The-QT-Bunny/tree-structrued-parser


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# @Time     : 2021/4/26-14:33
# @Author   : TuringEmmy
# @Email    : yonglonggeng@163.com
# @WeChat   : superior_god
# @File     : server.py
# @Project  : 00PythonProjects
import argparse
import os

import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf

from utils import montage

parser = argparse.ArgumentParser()
parser.add_argument('--auther', help='inner batch size', default="TuringEmmy", type=str)
parser.add_argument('--OUTPUT', help='weather to show graph', default=".", type=str, required=True)
parser.add_argument('--BATCH_SIZE', help='weather to show graph', default=128, type=int)
parser.add_argument('--Z_DIM', help='weather to show graph', default=100, type=int)
parser.add_argument('--DATASETS_NAME', help='lfw_new_imgs  or celeba datasets', default="lfw_new_imgs", type=str)

args = parser.parse_args()

with tf.Session() as sess:
    sess.run(tf.global_variables_initializer())

    model_path = args.OUTPUT + '/samples_' + args.DATASETS_NAME
    saver = tf.train.import_meta_graph(
        os.path.join(model_path, 'dcgan_' + args.DATASETS_NAME + '-60000.meta'))
    saver.restore(sess, tf.train.latest_checkpoint(model_path))
    graph = tf.get_default_graph()
    g = graph.get_tensor_by_name('generator/g/Tanh:0')
    noise = graph.get_tensor_by_name('noise:0')
    is_training = graph.get_tensor_by_name('is_training:0')

    n = np.random.uniform(-1.0, 1.0, [args.BATCH_SIZE, args.Z_DIM]).astype(np.float32)
    gen_imgs = sess.run(g, feed_dict={noise: n, is_training: False})
    gen_imgs = (gen_imgs + 1) / 2
    imgs = [img[:, :, :] for img in gen_imgs]
    gen_imgs = montage(imgs)
    gen_imgs = np.clip(gen_imgs, 0, 1)
    plt.figure(figsize=(8, 8))
    plt.axis('off')
    plt.imshow(gen_imgs)
    plt.show()

max_stars_repo_name: DragonYong/DCGANFACEGENERATION

max_issues_repo_name: DragonYong/DCGANFACEGENERATION

max_forks_repo_name: DragonYong/DCGANFACEGENERATION


text:
import numpy as np
import torch
import torch.nn.functional as F
from torch import nn

# https://github.com/adobe/antialiased-cnns
def get_pad_layer(pad_type):
    if (pad_type in ['refl', 'reflect']):
        PadLayer = nn.ReflectionPad2d
    elif (pad_type in ['repl', 'replicate']):
        PadLayer = nn.ReplicationPad2d
    elif (pad_type == 'zero'):
        PadLayer = nn.ZeroPad2d
    else:
        print('Pad type [%s] not recognized' % pad_type)
    return PadLayer

class Downsample(nn.Module):
    def __init__(self, pad_type='reflect', filt_size=3, stride=2, channels=None, pad_off=0):
        super(Downsample, self).__init__()
        self.filt_size = filt_size
        self.pad_off = pad_off
        self.pad_sizes = [int(1. * (filt_size - 1) / 2), int(np.ceil(1. * (filt_size - 1) / 2)),
                          int(1. * (filt_size - 1) / 2), int(np.ceil(1. * (filt_size - 1) / 2))]
        self.pad_sizes = [pad_size + pad_off for pad_size in self.pad_sizes]
        self.stride = stride
        self.off = int((self.stride - 1) / 2.)
        self.channels = channels

        # print('Filter size [%i]'%filt_size)
        assert self.filt_size > 0 and self.filt_size < 8
        a_list = [[], [1., ], [1., 1.], [1., 2., 1.], [1., 3., 3., 1.], [1., 4., 6., 4., 1.],
                  [1., 5., 10., 10., 5., 1.], [1., 6., 15., 20., 15., 6., 1.]]
        a = np.array(a_list[self.filt_size])

        filt = torch.Tensor(a[:, None] * a[None, :])
        filt = filt / torch.sum(filt)
        self.register_buffer('filt', filt[None, None, :, :].repeat((self.channels, 1, 1, 1)))

        self.pad = get_pad_layer(pad_type)(self.pad_sizes)

    def forward(self, inp):
        if (self.filt_size == 1):
            if (self.pad_off == 0):
                return inp[:, :, ::self.stride, ::self.stride]
            else:
                return self.pad(inp)[:, :, ::self.stride, ::self.stride]
        else:
            return F.conv2d(self.pad(inp), self.filt, stride=self.stride, groups=inp.shape[1])

class MaxBlurPool2d(nn.Module):
    def __init__(self, kernel_size, stride=None, padding=0, dilation=1, return_indices=False, ceil_mode=False,
                 channels=None, filt_size=5):
        super(MaxBlurPool2d, self).__init__()
        self.kernel_size = kernel_size
        self.stride = stride or kernel_size
        self.padding = padding
        self.dilation = dilation
        self.return_indices = return_indices
        self.ceil_mode = ceil_mode

        if stride == 1:
            self.pool = nn.MaxPool2d(kernel_size, stride, padding, dilation, return_indices, ceil_mode)
        else:
            self.pool = nn.Sequential(nn.MaxPool2d(kernel_size, 1, padding, dilation, return_indices, ceil_mode),
                                      Downsample(channels=channels, filt_size=filt_size, stride=stride))

    def forward(self, x):
        return self.pool(x)


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import os
import sys
import random
import numpy as np
import tensorflow as tf
from PIL import Image
from datasets import utils

# The height and width of each image.
_IMAGE_SIZE = 32

def _write_to_tfrecord(filenames, labels, tfrecord_writer):
    num_images = len(filenames)
    with tf.Graph().as_default():
        image_placeholder = tf.placeholder(dtype=tf.uint8)
        encoded_image = tf.image.encode_png(image_placeholder)
        with tf.Session('') as sess:
            for i in range(num_images):
                sys.stdout.write('\r>> Reading images %d/%d' % (i + 1, num_images))
                sys.stdout.flush()
                image_path = filenames[i]
                image = Image.open(image_path)
                label = labels[i]
                png_string = sess.run(encoded_image, feed_dict={image_placeholder: image})
                example = utils.image_to_tfexample(png_string, 'png', _IMAGE_SIZE, _IMAGE_SIZE, label)
                tfrecord_writer.write(example.SerializeToString())

def _get_testdata_filenames_and_labels(dataset_dir):
    classes = os.listdir(dataset_dir)
    filenames = []
    labels = []
    for i in range(len(classes)):
        classpath = os.path.join(os.getcwd(), dataset_dir, classes[i])
        imgpaths = os.listdir(classpath)
        imgpaths = [os.path.join(classpath, substring) for substring in imgpaths]
        filenames.append(imgpaths)
        labels.append(i * np.ones(len(imgpaths), dtype=np.int64))
    filenames = np.concatenate(filenames, 0)
    labels = np.concatenate(labels, 0)
    return filenames, labels

def _shuffle_data(filenames, labels):
    shuffled_index = list(range(len(filenames)))
    random.seed(12345)
    random.shuffle(shuffled_index)
    filenames = [filenames[i] for i in shuffled_index]
    labels = [labels[i] for i in shuffled_index]
    return filenames, labels

def _get_traindata_filenames_and_labels(dataset_dir, labels_per_class):
    classes = os.listdir(dataset_dir)
    filenames_l = []
    filenames_u = []
    labels_l = []
    labels_u = []
    for i in range(len(classes)):
        classpath = os.path.join(os.getcwd(), dataset_dir, classes[i])
        imgpaths = os.listdir(classpath)
        imgpaths = [os.path.join(classpath, substring) for substring in imgpaths]
        random.seed(12345)
        random.shuffle(imgpaths)
        imgpaths_l = imgpaths[:labels_per_class]
        imgpaths_u = imgpaths[labels_per_class:]
        filenames_l.append(imgpaths_l)
        filenames_u.append(imgpaths_u)
        labels_l.append(i * np.ones(labels_per_class, dtype=np.int64))
        labels_u.append(-1 * np.ones(len(imgpaths) - labels_per_class, dtype=np.int64))
    filenames_l = np.concatenate(filenames_l, 0)
    filenames_u = np.concatenate(filenames_u, 0)
    labels_l = np.concatenate(labels_l, 0)
    labels_u = np.concatenate(labels_u, 0)
    # shuffle traindata completely
    filenames_l, labels_l = _shuffle_data(filenames_l, labels_l)
    filenames_u, labels_u = _shuffle_data(filenames_u, labels_u)
    return filenames_l, labels_l, filenames_u, labels_u

def _make_dir(output_dir):
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

def _get_tfrecords_names(output_dir, dataset_name, split_name):
    tfrecords_names = '%s/%s_%s.tfrecord' % (output_dir, dataset_name, split_name)
    _make_dir(output_dir)
    return tfrecords_names

def run(dataset_dir, split_name, dataset_name, output_dir, num_splits, labels_per_class):

    if split_name == 'test':
        tfrecords_names = _get_tfrecords_names(output_dir, dataset_name, split_name)
        filenames, labels = _get_testdata_filenames_and_labels(dataset_dir)
        with tf.python_io.TFRecordWriter(tfrecords_names) as tfrecord_writer:
            _write_to_tfrecord(filenames, labels, tfrecord_writer)
    else:
        for i in range(num_splits):
            output_sub_dir_l = os.path.join(output_dir, 'split'+str(i+1), 'labeled')
            tfrecords_names_l = _get_tfrecords_names(output_sub_dir_l, dataset_name, split_name)

            output_sub_dir_u = os.path.join(output_dir, 'split'+str(i+1), 'unlabeled')
            tfrecords_names_u = _get_tfrecords_names(output_sub_dir_u, dataset_name, split_name)

            filenames_l, labels_l, filenames_u, labels_u = \
                _get_traindata_filenames_and_labels(dataset_dir, labels_per_class)

            print('\nConverting the labelled data.')
            with tf.python_io.TFRecordWriter(tfrecords_names_l) as tfrecord_writer:
                _write_to_tfrecord(filenames_l, labels_l, tfrecord_writer)

            print('\nConverting the unlabelled data.')
            with tf.python_io.TFRecordWriter(tfrecords_names_u) as tfrecord_writer:
                _write_to_tfrecord(filenames_u, labels_u, tfrecord_writer)

    print('\nFinished converting the data.')

max_stars_repo_name: yanbeic/semi-memory

max_issues_repo_name: Matanley/semi-memory

max_forks_repo_name: Matanley/semi-memory


text:
# by Huaiyu-TAN

# m: initial number of nodes
# n: final number of nodes
# k: nodes degree

import networkx as nx
import numpy as np
from time import time

class BA_sfnet:

    def __init__(self, m, n, k):
        self.m = m; self.n = n; self.k = k
        self.G = nx.complete_graph(self.m)
        self.new_node = self.m
    def choice_node(self):
        # print("最初节点有 {} 个".format(len(self.G.nodes)))
        nodes_probs = []
        for node in self.G.nodes():
            nodes_probs.append(self.G.degree(node) / (2 * len(self.G.edges()))) # 连边概率
        # print("连边概率：{}".format(nodes_probs))
        random_prob_node = np.random.choice(self.G.nodes(), p = nodes_probs) # 以连边概率p从全部节点中选择一个与新节点相连
        # print("随机区间内选择节点 {} 号".format(random_prob_node))
        return random_prob_node

    def add_edge(self): # 采用递归
        if len(self.G.edges()) == 0: # 基线条件
            random_prob_node = 0
        else:
            random_prob_node = BA_sfnet.choice_node(self)
        new_edge = (random_prob_node, self.new_node)
        if new_edge in self.G.edges():
            self.G.add_edge
        else:
            self.G.add_edge(self.new_node, random_prob_node) # 加入新边
            # print("加入边: ({} {})".format(self.new_node, random_prob_node)) # new_node+1显示错位

    def sfnet(self):
        count = 0
        start = time()
        for i in range(self.n - self.m):
            print("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<第 {} 轮模拟>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>".format(count + 1))
            self.G.add_node(self.m + count)
            # print("添加第 {} 个节点".format(self.m + count + 1))
            count += 1
            for j in range(0, self.k):
                BA_sfnet.add_edge(self)
            self.new_node += 1
        end = time()

        # python排序从0开始，这里统一加1
        list_edges = list(self.G.edges)
        # edges_point1 = [list(G.edges)[i][0] + 1 for i in range(len(list_edges))]
        # edges_point2 = [list(G.edges)[i][1] + 1 for i in range(len(list_edges))]

        # 储存为邻接矩阵
        nebmatrix = np.zeros(shape = (self.n, self.n))
        for i in range(len(list_edges)):
            nebmatrix[ list_edges[i][0], list_edges[i][1] ] = 1
            nebmatrix[ list_edges[i][1], list_edges[i][0] ] = 1
        # np.savetxt('Adjacent_Matrix.csv', nebmatrix, delimiter = ',')
        print(nebmatrix)
        # print("\n\n 最终{}个节点的连边情况为：\n\n {} \n".format(self.n, self.G.edges()))
        # print("\n 系统内存在的节点：\n\n{}\n".format(self.G.nodes()))
        print("\n \n {}轮用时共为{}秒".format(self.n ,end - start))

max_stars_repo_path: Some Complex Network simulation/BA_scale-free.py

reason: import numpy,import networkx


text:
import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F
import dgl
from dgl.nn import GraphConv
from models.AED.simpleAED import Encoder
from models.discriminator import Discriminator
from models.losses import global_loss_

class GraphEncoder(nn.Module):
    def __init__(self, node_nums, in_feats, h1_feats, last_space_feature):
        super(GraphEncoder, self).__init__()
        # self.conv1 = GATConv(in_feats, h1_feats, 1)
        self.node_nums = node_nums
        self.in_feats = in_feats
        self.last_space_feature = last_space_feature
        self.conv1 = GraphConv(in_feats, h1_feats)
        self.conv2 = GraphConv(h1_feats, last_space_feature)
        # self.conv2 = GATConv(h1_feats, last_space_feature, 1)
        # self.device = device
        self.init_emb()

    def init_emb(self):
        for m in self.modules():
            if isinstance(m, nn.Linear):
                torch.nn.init.xavier_uniform_(m.weight.data)
                if m.bias is not None:
                    m.bias.data.fill_(0.0)

    def forward(self, g):
        batch_size = int(g.num_nodes() / self.node_nums)
        h = self.conv1(g, g.ndata['features'])
        h = F.leaky_relu(h)
        h1 = h
        g.ndata['h1'] = h1
        h = self.conv2(g, h)
        h = F.leaky_relu(h)
        h2 = h
        # H = torch.reshape(h2, (batch_size, self.node_nums, self.last_space_feature))
        # H = torch.mean(H, axis=1)
        g.ndata['h2'] = h2
        H = dgl.readout.mean_nodes(g, 'h2')
        # g.ndata['H'] = H
        return H

class GraphPMU(nn.Module):
    def __init__(self, encoder, discriminator, node_nums, in_feats, h1_feats, last_space_feature, D_h1, D_h2, device):
        super(GraphPMU, self).__init__()
        self.node_nums = node_nums
        self.in_feats = in_feats
        self.h1_feats = h1_feats
        self.last_space_feature = last_space_feature
        self.D_h1 = D_h1
        self.D_h2 = D_h2
        # self.measure = measure
        self.device = device
        self.encoder = encoder(node_nums, in_feats, h1_feats, last_space_feature).to(device)
        # self.discriminator = discriminator(last_space_feature, D_h1, D_h2).to(device)#if simple encoder
        self.discriminator = discriminator(last_space_feature + h1_feats, D_h1, D_h2).to(device)#if locglob encoder
        self.init_emb()

    def init_emb(self):
        for m in self.modules():
            if isinstance(m, nn.Linear):
                torch.nn.init.xavier_uniform_(m.weight.data)
                if m.bias is not None:
                    m.bias.data.fill_(0.0)

    def forward(self, g):
        H = self.encoder(g.to(self.device))
        y = self.discriminator(H)

        return y

class GraphEncoderLocGlob(nn.Module):
    def __init__(self, node_nums, in_feats, h1_feats, last_space_feature):
        super(GraphEncoderLocGlob, self).__init__()
        # self.conv1 = GATConv(in_feats, h1_feats, 1)
        self.node_nums = node_nums
        self.in_feats = in_feats
        self.last_space_feature = last_space_feature
        self.conv1 = GraphConv(in_feats, h1_feats)
        self.conv2 = GraphConv(h1_feats, last_space_feature)
        # self.conv2 = GATConv(h1_feats, last_space_feature, 1)
        # self.device = device
        self.init_emb()

    def init_emb(self):
        for m in self.modules():
            if isinstance(m, nn.Linear):
                torch.nn.init.xavier_uniform_(m.weight.data)
                if m.bias is not None:
                    m.bias.data.fill_(0.0)

    def forward(self, g):
        batch_size = int(g.num_nodes() / self.node_nums)
        h = self.conv1(g, g.ndata['features'])
        h = F.leaky_relu(h)
        h1 = h
        g.ndata['h1'] = h1
        h = self.conv2(g, h)
        h = F.leaky_relu(h)
        h2 = h
        # H = torch.reshape(h2, (batch_size, self.node_nums, self.last_space_feature))
        # H = torch.mean(H, axis=1)
        g.ndata['h2'] = h2
        g.ndata['hcat'] = torch.cat((h1, h2), 1)
        H = dgl.readout.mean_nodes(g, 'hcat')
        # g.ndata['H'] = H
        return H

class GraphPMULocalGlobal(nn.Module):
    def __init__(self, encoder, discriminator, node_nums, in_feats, h1_feats, last_space_feature, D_h1, D_h2, device):
        super(GraphPMULocalGlobal, self).__init__()
        self.node_nums = node_nums
        self.in_feats = in_feats
        self.h1_feats = h1_feats
        self.last_space_feature = last_space_feature
        self.D_h1 = D_h1
        self.D_h2 = D_h2
        # self.measure = measure
        self.device = device
        self.encoder = encoder(node_nums, in_feats, h1_feats, last_space_feature).to(device)
        # self.discriminator = discriminator(last_space_feature, D_h1, D_h2).to(device)#if simple encoder
        # self.discriminator = discriminator(2 * (last_space_feature + h1_feats), D_h1, D_h2).to(device)#if locglob encoder
        self.discriminator = discriminator(last_space_feature + h1_feats, D_h1, D_h2).to(device)#if locglob encoder
        self.init_emb()

    def init_emb(self):
        for m in self.modules():
            if isinstance(m, nn.Linear):
                torch.nn.init.xavier_uniform_(m.weight.data)
                if m.bias is not None:
                    m.bias.data.fill_(0.0)

    def forward(self, g):
        H = self.encoder(g.to(self.device)) # encoder should be GraphEncoderLocGlob
        y = self.discriminator(H)

        return y

class GraphEncoderLocGlobAutoEncoder(nn.Module):
    def __init__(self, node_nums, in_feats, h1_feats, last_space_feature):
        super(GraphEncoderLocGlobAutoEncoder, self).__init__()
        # self.conv1 = GATConv(in_feats, h1_feats, 1)
        self.node_nums = node_nums
        self.in_feats = in_feats
        self.last_space_feature = last_space_feature
        self.conv1 = GraphConv(in_feats, h1_feats)
        self.conv2 = GraphConv(h1_feats, last_space_feature)
        # self.conv2 = GATConv(h1_feats, last_space_feature, 1)
        # self.device = device
        self.init_emb()

    def init_emb(self):
        for m in self.modules():
            if isinstance(m, nn.Linear):
                torch.nn.init.xavier_uniform_(m.weight.data)
                if m.bias is not None:
                    m.bias.data.fill_(0.0)

    def forward(self, g):
        batch_size = int(g.num_nodes() / self.node_nums)
        h = self.conv1(g, g.ndata['embd'])
        h = F.leaky_relu(h)
        h1 = h
        g.ndata['h1'] = h1
        h = self.conv2(g, h)
        h = F.leaky_relu(h)
        h2 = h
        # H = torch.reshape(h2, (batch_size, self.node_nums, self.last_space_feature))
        # H = torch.mean(H, axis=1)
        g.ndata['h2'] = h2
        g.ndata['hcat'] = torch.cat((h1, h2), 1)
        H = dgl.readout.mean_nodes(g, 'hcat')
        # g.ndata['H'] = H
        return g.ndata['hcat'], H

class EncGraphDisc(nn.Module):
    def __init__(self, encoder,graphencoder, discriminator, seq_len, n_time_features,enc_emb_size, node_nums, h1_feats,
                 last_space_feature, D_h1, D_h2, device):
        super(EncGraphDisc, self).__init__()

        self.seq_len = seq_len
        self.n_time_features = n_time_features
        self.node_nums = node_nums
        self.enc_emb_size = enc_emb_size
        self.h1_feats = h1_feats
        self.last_space_feature = last_space_feature
        self.D_h1 = D_h1
        self.D_h2 = D_h2
        # self.measure = measure
        self.device = device
        self.encoder = encoder(seq_len, n_time_features, enc_emb_size).to(device)
        self.graphencoder = graphencoder(node_nums, enc_emb_size, h1_feats
                                                           , last_space_feature).to(device)
        # self.discriminator = discriminator(last_space_feature, D_h1, D_h2).to(device)#if simple encoder
        self.discriminator = discriminator(2 * (last_space_feature + h1_feats), D_h1, D_h2).to(device)#if locglob encoder
        self.init_emb()

    def init_emb(self):
        for m in self.modules():
            if isinstance(m, nn.Linear):
                torch.nn.init.xavier_uniform_(m.weight.data)
                if m.bias is not None:
                    m.bias.data.fill_(0.0)

    def forward(self, g):
        g = g.to(self.device)
        g.ndata['embd'] = self.encoder(g.ndata['features'].to(self.device))
        g.ndata['hcat'], H = self.graphencoder(g.to(self.device))  # encoder should be GraphEncoderLocGlobAutoEncoder
        # y = self.discriminator(H)
        return g.ndata['embd'], g.ndata['hcat'], H

max_stars_repo_name: arminalgln/graphpmu

max_issues_repo_name: arminalgln/graphpmu

max_forks_repo_name: arminalgln/graphpmu


text:
import numpy as np
import pandas as pd

a = np.random.randn(2,3)
print(a)

b = np.random.random((2,3))
print(b)

max_stars_repo_path: 01-python/source code/03/07.py

max_issues_repo_path: 01-python/source code/03/07.py

max_forks_repo_path: 01-python/source code/03/07.py