text:
#!/usr/bin/python
#
# Copyright 2022 DeepMind Technologies Limited
#
# 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.

"""Common Losses to be used to train a model."""

import chex
import jax.numpy as jnp
import optax

def softmax_cross_entropy(logits: chex.Array,
                          labels: chex.Array,
                          reduction: str = 'sum') -> chex.Array:
  """Computes softmax cross entropy given logits and one-hot class labels.

  Args:
    logits: Logit output values.
    labels: Ground truth one-hot-encoded labels.
    reduction: Type of reduction to apply to loss.

  Returns:
    Loss value. If `reduction` is `none`, this has the same shape as `labels`;
    otherwise, it is scalar.

  Raises:
    ValueError: If the type of `reduction` is unsupported.
  """
  x = optax.softmax_cross_entropy(logits, labels)
  if reduction == 'none' or reduction is None:
    return jnp.asarray(x)
  elif reduction == 'sum':
    return jnp.asarray(x).sum()
  elif reduction == 'mean':
    return jnp.mean(jnp.asarray(x))
  else:
    raise ValueError('Unsupported reduction option.')

reason: import jax


text:
#-------------------------------------------------------------
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you 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.
#
#-------------------------------------------------------------

__all__ = [ 'setSparkContext', 'matrix', 'eval', 'solve', 'DMLOp', 'set_lazy', 'debug_array_conversion', 'load', 'full', 'seq' ]

import numpy as np
import pandas as pd
from scipy.sparse import coo_matrix, spmatrix
try:
    import py4j.java_gateway
    from py4j.java_gateway import JavaObject
    from pyspark import SparkContext
    from pyspark.sql import DataFrame, SQLContext
    import pyspark.mllib.common
except ImportError:
    raise ImportError('Unable to import `pyspark`. Hint: Make sure you are running with PySpark.')

from . import MLContext, pydml, _java2py, Matrix
from .converters import *

def setSparkContext(sc):
    """
    Before using the matrix, the user needs to invoke this function if SparkContext is not previously created in the session.

    Parameters
    ----------
    sc: SparkContext
        SparkContext
    """
    matrix.sc = sc
    matrix.sqlContext = SQLContext(sc)
    matrix.ml = MLContext(matrix.sc)

def check_MLContext():
    if matrix.ml is None:
        if SparkContext._active_spark_context is not None:
            setSparkContext(SparkContext._active_spark_context)
        else:
            raise Exception('Expected setSparkContext(sc) to be called, where sc is active SparkContext.')

########################## AST related operations ##################################

class DMLOp(object):
    """
    Represents an intermediate node of Abstract syntax tree created to generate the PyDML script
    """
    def __init__(self, inputs, dml=None):
        self.inputs = inputs
        self.dml = dml
        self.ID = None
        self.depth = 1
        for m in self.inputs:
            m.referenced = m.referenced + [ self ]
            if isinstance(m, matrix) and m.op is not None:
                self.depth = max(self.depth, m.op.depth + 1)

    MAX_DEPTH = 0

    def _visit(self, execute=True):
        matrix.dml = matrix.dml + self.dml

    def _print_ast(self, numSpaces):
        ret = []
        for m in self.inputs:
            ret = [ m._print_ast(numSpaces+2) ]
        return ''.join(ret)

# Special object used internally to specify the placeholder which will be replaced by output ID
# This helps to provide dml containing output ID in construct_intermediate_node
OUTPUT_ID = '$$OutputID$$'

def set_lazy(isLazy):
    """
    This method allows users to set whether the matrix operations should be executed in lazy manner.

    Parameters
    ----------
    isLazy: True if matrix operations should be evaluated in lazy manner.
    """
    if isLazy:
        DMLOp.MAX_DEPTH = 0
    else:
        DMLOp.MAX_DEPTH = 1

def construct_intermediate_node(inputs, dml):
    """
    Convenient utility to create an intermediate node of AST.

    Parameters
    ----------
    inputs = list of input matrix objects and/or DMLOp
    dml = list of DML string (which will be eventually joined before execution). To specify out.ID, please use the placeholder
    """
    dmlOp = DMLOp(inputs)
    out = matrix(None, op=dmlOp)
    dmlOp.dml = [out.ID if x==OUTPUT_ID else x for x in dml]
    if DMLOp.MAX_DEPTH > 0 and out.op.depth >= DMLOp.MAX_DEPTH:
        out.eval()
    return out

def load(file, format='csv'):
    """
    Allows user to load a matrix from filesystem

    Parameters
    ----------
    file: filepath
    format: can be csv, text or binary or mm
    """
    return construct_intermediate_node([], [OUTPUT_ID, ' = load(\"', file, '\", format=\"', format, '\")\n'])

def full(shape, fill_value):
    """
    Return a new array of given shape filled with fill_value.

    Parameters
    ----------
    shape: tuple of length 2
    fill_value: float or int
    """
    return construct_intermediate_node([], [OUTPUT_ID, ' = full(', str(fill_value), ', rows=', str(shape[0]), ', cols=', str(shape[1]), ')\n'])

def reset():
    """
    Resets the visited status of matrix and the operators in the generated AST.
    """
    for m in matrix.visited:
        m.visited = False
    matrix.visited = []
    matrix.ml = MLContext(matrix.sc)
    matrix.dml = []
    matrix.script = pydml('')

def perform_dfs(outputs, execute):
    """
    Traverses the forest of nodes rooted at outputs nodes and returns the DML script to execute
    """
    for m in outputs:
        m.output = True
        m._visit(execute=execute)
    return ''.join(matrix.dml)

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

########################## Utility functions ##################################

def _log_base(val, base):
    if not isinstance(val, str):
        raise ValueError('The val to _log_base should be of type string')
    return '(log(' + val + ')/log(' + str(base) + '))'

def _matricize(lhs, inputs):
    """
    Utility fn to convert the supported types to matrix class or to string (if float or int)
    and return the string to be passed to DML as well as inputs
    """
    if isinstance(lhs, SUPPORTED_TYPES):
        lhs = matrix(lhs)
    if isinstance(lhs, matrix):
        lhsStr = lhs.ID
        inputs = inputs + [lhs]
    elif isinstance(lhs, float) or isinstance(lhs, int):
        lhsStr = str(lhs)
    else:
        raise TypeError('Incorrect type')
    return lhsStr, inputs

def binary_op(lhs, rhs, opStr):
    """
    Common function called by all the binary operators in matrix class
    """
    inputs = []
    lhsStr, inputs = _matricize(lhs, inputs)
    rhsStr, inputs = _matricize(rhs, inputs)
    return construct_intermediate_node(inputs, [OUTPUT_ID, ' = ', lhsStr, opStr, rhsStr, '\n'])

def binaryMatrixFunction(X, Y, fnName):
    """
    Common function called by supported PyDML built-in function that has two arguments.
    """
    inputs = []
    lhsStr, inputs = _matricize(X, inputs)
    rhsStr, inputs = _matricize(Y, inputs)
    return construct_intermediate_node(inputs, [OUTPUT_ID, ' = ', fnName,'(', lhsStr, ', ', rhsStr, ')\n'])

def unaryMatrixFunction(X, fnName):
    """
    Common function called by supported PyDML built-in function that has one argument.
    """
    inputs = []
    lhsStr, inputs = _matricize(X, inputs)
    return construct_intermediate_node(inputs, [OUTPUT_ID, ' = ', fnName,'(', lhsStr, ')\n'])

def seq(start=None, stop=None, step=1):
    """
    Creates a single column vector with values starting from <start>, to <stop>, in increments of <step>.
    Note: Unlike Numpy's arange which returns a row-vector, this returns a column vector.
    Also, Unlike Numpy's arange which doesnot include stop, this method includes stop in the interval.

    Parameters
    ----------
    start: int or float [Optional: default = 0]
    stop: int or float
    step : int float [Optional: default = 1]
    """
    if start is None and stop is None:
        raise ValueError('Both start and stop cannot be None')
    elif start is not None and stop is None:
        stop = start
        start = 0
    return construct_intermediate_node([], [OUTPUT_ID, ' = seq(', str(start), ',', str(stop), ',',  str(step), ')\n'])

# utility function that converts 1:3 into DML string
def convert_seq_to_dml(s):
    ret = []
    if s is None:
        return ''
    elif isinstance(s, slice):
        if s.step is not None:
            raise ValueError('Slicing with step is not supported.')
        if s.start is None:
            ret = ret + [ '0 : ' ]
        else:
            ret = ret + [ getValue(s.start), ':' ]
        if s.start is None:
            ret = ret + [ '' ]
        else:
            ret = ret + [ getValue(s.stop) ]
    else:
        ret = ret + [ getValue(s) ]
    return ''.join(ret)

# utility function that converts index (such as [1, 2:3]) into DML string
def getIndexingDML(index):
    ret = [ '[' ]
    if isinstance(index, tuple) and len(index) == 1:
        ret = ret + [ convert_seq_to_dml(index[0]), ',' ]
    elif isinstance(index, tuple) and len(index) == 2:
        ret = ret + [ convert_seq_to_dml(index[0]), ',', convert_seq_to_dml(index[1]) ]
    else:
        raise TypeError('matrix indexes can only be tuple of length 2. For example: m[1,1], m[0:1,], m[:, 0:1]')
    return ret + [ ']' ]

def convert_outputs_to_list(outputs):
    if isinstance(outputs, matrix):
        return [ outputs ]
    elif isinstance(outputs, list):
        for o in outputs:
            if not isinstance(o, matrix):
                raise TypeError('Only matrix or list of matrix allowed')
        return outputs
    else:
        raise TypeError('Only matrix or list of matrix allowed')

def reset_output_flag(outputs):
    for m in outputs:
        m.output = False

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

########################## Global user-facing functions #######################

def solve(A, b):
    """
    Computes the least squares solution for system of linear equations A %*% x = b

    Examples
    --------
    >>> import numpy as np
    >>> from sklearn import datasets
    >>> import SystemML as sml
    >>> from pyspark.sql import SQLContext
    >>> diabetes = datasets.load_diabetes()
    >>> diabetes_X = diabetes.data[:, np.newaxis, 2]
    >>> X_train = diabetes_X[:-20]
    >>> X_test = diabetes_X[-20:]
    >>> y_train = diabetes.target[:-20]
    >>> y_test = diabetes.target[-20:]
    >>> sml.setSparkContext(sc)
    >>> X = sml.matrix(X_train)
    >>> y = sml.matrix(y_train)
    >>> A = X.transpose().dot(X)
    >>> b = X.transpose().dot(y)
    >>> beta = sml.solve(A, b).toNumPy()
    >>> y_predicted = X_test.dot(beta)
    >>> print('Residual sum of squares: %.2f' % np.mean((y_predicted - y_test) ** 2))
    Residual sum of squares: 25282.12
    """
    return binaryMatrixFunction(A, b, 'solve')

def eval(outputs, execute=True):
    """
    Executes the unevaluated DML script and computes the matrices specified by outputs.

    Parameters
    ----------
    outputs: list of matrices or a matrix object
    execute: specified whether to execute the unevaluated operation or just return the script.
    """
    check_MLContext()
    reset()
    outputs = convert_outputs_to_list(outputs)
    matrix.script.scriptString = perform_dfs(outputs, execute)
    if not execute:
        reset_output_flag(outputs)
        return matrix.script.scriptString
    results = matrix.ml.execute(matrix.script)
    for m in outputs:
        m.eval_data = results._java_results.get(m.ID)
    reset_output_flag(outputs)

def debug_array_conversion(throwError):
    matrix.THROW_ARRAY_CONVERSION_ERROR = throwError

def _get_new_var_id():
    matrix.systemmlVarID += 1
    return 'mVar' + str(matrix.systemmlVarID)

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

class matrix(object):
    """
    matrix class is a python wrapper that implements basic matrix operators, matrix functions
    as well as converters to common Python types (for example: Numpy arrays, PySpark DataFrame
    and Pandas DataFrame).

    The operators supported are:

    1. Arithmetic operators: +, -, *, /, //, %, ** as well as dot (i.e. matrix multiplication)
    2. Indexing in the matrix
    3. Relational/Boolean operators: <, <=, >, >=, ==, !=, &, |

    In addition, following functions are supported for matrix:

    1. transpose
    2. Aggregation functions: sum, mean, var, sd, max, min, argmin, argmax, cumsum
    3. Global statistical built-In functions: exp, log, abs, sqrt, round, floor, ceil, sin, cos, tan, asin, acos, atan, sign, solve

    For all the above functions, we always return a two dimensional matrix, especially for aggregation functions with axis.
    For example: Assuming m1 is a matrix of (3, n), NumPy returns a 1d vector of dimension (3,) for operation m1.sum(axis=1)
    whereas SystemML returns a 2d matrix of dimension (3, 1).

    Note: an evaluated matrix contains a data field computed by eval method as DataFrame or NumPy array.

    Examples
    --------
    >>> import SystemML as sml
    >>> import numpy as np
    >>> sml.setSparkContext(sc)

    Welcome to Apache SystemML!

    >>> m1 = sml.matrix(np.ones((3,3)) + 2)
    >>> m2 = sml.matrix(np.ones((3,3)) + 3)
    >>> m2 = m1 * (m2 + m1)
    >>> m4 = 1.0 - m2
    >>> m4
    # This matrix (mVar5) is backed by below given PyDML script (which is not yet evaluated). To fetch the data of this matrix, invoke toNumPy() or toDF() or toPandas() methods.
    mVar1 = load(" ", format="csv")
    mVar2 = load(" ", format="csv")
    mVar3 = mVar2 + mVar1
    mVar4 = mVar1 * mVar3
    mVar5 = 1.0 - mVar4
    save(mVar5, " ")
    >>> m2.eval()
    >>> m2
    # This matrix (mVar4) is backed by NumPy array. To fetch the NumPy array, invoke toNumPy() method.
    >>> m4
    # This matrix (mVar5) is backed by below given PyDML script (which is not yet evaluated). To fetch the data of this matrix, invoke toNumPy() or toDF() or toPandas() methods.
    mVar4 = load(" ", format="csv")
    mVar5 = 1.0 - mVar4
    save(mVar5, " ")
    >>> m4.sum(axis=1).toNumPy()
    array([[-60.],
           [-60.],
           [-60.]])

    Design Decisions:

    1. Until eval() method is invoked, we create an AST (not exposed to the user) that consist of unevaluated operations and data required by those operations.
       As an anology, a spark user can treat eval() method similar to calling RDD.persist() followed by RDD.count().
    2. The AST consist of two kinds of nodes: either of type matrix or of type DMLOp.
       Both these classes expose _visit method, that helps in traversing the AST in DFS manner.
    3. A matrix object can either be evaluated or not.
       If evaluated, the attribute 'data' is set to one of the supported types (for example: NumPy array or DataFrame). In this case, the attribute 'op' is set to None.
       If not evaluated, the attribute 'op' which refers to one of the intermediate node of AST and if of type DMLOp.  In this case, the attribute 'data' is set to None.
    5. DMLOp has an attribute 'inputs' which contains list of matrix objects or DMLOp.
    6. To simplify the traversal, every matrix object is considered immutable and an matrix operations creates a new matrix object.
       As an example:
       `m1 = sml.matrix(np.ones((3,3)))` creates a matrix object backed by 'data=(np.ones((3,3))'.
       `m1 = m1 * 2` will create a new matrix object which is now backed by 'op=DMLOp( ... )' whose input is earlier created matrix object.
    7. Left indexing (implemented in __setitem__ method) is a special case, where Python expects the existing object to be mutated.
       To ensure the above property, we make deep copy of existing object and point any references to the left-indexed matrix to the newly created object.
       Then the left-indexed matrix is set to be backed by DMLOp consisting of following pydml:
       left-indexed-matrix = new-deep-copied-matrix
       left-indexed-matrix[index] = value
    8. Please use m.print_ast() and/or  type `m` for debugging. Here is a sample session:

       >>> npm = np.ones((3,3))
       >>> m1 = sml.matrix(npm + 3)
       >>> m2 = sml.matrix(npm + 5)
       >>> m3 = m1 + m2
       >>> m3
       mVar2 = load(" ", format="csv")
       mVar1 = load(" ", format="csv")
       mVar3 = mVar1 + mVar2
       save(mVar3, " ")
       >>> m3.print_ast()
       - [mVar3] (op).
         - [mVar1] (data).
         - [mVar2] (data).

    """
    # Global variable that is used to keep track of intermediate matrix variables in the DML script
    systemmlVarID = 0

    # Since joining of string is expensive operation, we collect the set of strings into list and then join
    # them before execution: See matrix.script.scriptString = ''.join(matrix.dml) in eval() method
    dml = []

    # Represents MLContext's script object
    script = None

    # Represents MLContext object
    ml = None

    # Contains list of nodes visited in Abstract Syntax Tree. This helps to avoid computation of matrix objects
    # that have been previously evaluated.
    visited = []

    def __init__(self, data, op=None):
        """
        Constructs a lazy matrix

        Parameters
        ----------
        data: NumPy ndarray, Pandas DataFrame, scipy sparse matrix or PySpark DataFrame. (data cannot be None for external users, 'data=None' is used internally for lazy evaluation).
        """
        self.dtype = np.double
        check_MLContext()
        self.visited = False
        self.output = False
        self.ID = _get_new_var_id()
        self.referenced = []
        # op refers to the node of Abstract Syntax Tree created internally for lazy evaluation
        self.op = op
        self.eval_data = data
        self._shape = None
        if isinstance(data, SUPPORTED_TYPES):
            self._shape = data.shape
        if not (isinstance(data, SUPPORTED_TYPES) or hasattr(data, '_jdf') or (data is None and op is not None)):
            raise TypeError('Unsupported input type')

    def eval(self):
        """
        This is a convenience function that calls the global eval method
        """
        eval([self])

    def toPandas(self):
        """
        This is a convenience function that calls the global eval method and then converts the matrix object into Pandas DataFrame.
        """
        self.eval()
        if isinstance(self.eval_data, py4j.java_gateway.JavaObject):
            self.eval_data = _java2py(SparkContext._active_spark_context, self.eval_data)
        if isinstance(self.eval_data, Matrix):
            self.eval_data = self.eval_data.toNumPy()
        self.eval_data = convertToPandasDF(self.eval_data)
        return self.eval_data

    def toNumPy(self):
        """
        This is a convenience function that calls the global eval method and then converts the matrix object into NumPy array.
        """
        self.eval()
        if isinstance(self.eval_data, py4j.java_gateway.JavaObject):
            self.eval_data = _java2py(SparkContext._active_spark_context, self.eval_data)
        if isinstance(self.eval_data, Matrix):
            self.eval_data = self.eval_data.toNumPy()
            return self.eval_data
        if isinstance(self.eval_data, pd.DataFrame):
            self.eval_data = self.eval_data.as_matrix()
        elif isinstance(self.eval_data, DataFrame):
            self.eval_data = self.eval_data.toPandas().as_matrix()
        elif isinstance(self.eval_data, spmatrix):
            self.eval_data = self.eval_data.toarray()
        elif isinstance(self.eval_data, Matrix):
            self.eval_data = self.eval_data.toNumPy()
        # Always keep default format as NumPy array if possible
        return self.eval_data

    def toDF(self):
        """
        This is a convenience function that calls the global eval method and then converts the matrix object into DataFrame.
        """
        if isinstance(self.eval_data, DataFrame):
            return self.eval_data
        if isinstance(self.eval_data, py4j.java_gateway.JavaObject):
            self.eval_data = _java2py(SparkContext._active_spark_context, self.eval_data)
        if isinstance(self.eval_data, Matrix):
            self.eval_data = self.eval_data.toDF()
            return self.eval_data
        self.eval_data = matrix.sqlContext.createDataFrame(self.toPandas())
        return self.eval_data

    def save(self, file, format='csv'):
        """
        Allows user to save a matrix to filesystem

        Parameters
        ----------
        file: filepath
        format: can be csv, text or binary or mm
        """
        tmp = construct_intermediate_node([self], ['save(', self.ID , ',\"', file, '\", format=\"', format, '\")\n'])
        construct_intermediate_node([tmp], [OUTPUT_ID, ' = full(0, rows=1, cols=1)\n']).eval()

    def _mark_as_visited(self):
        self.visited = True
        # for cleanup
        matrix.visited = matrix.visited + [ self ]
        return self

    def _register_as_input(self, execute):
        # TODO: Remove this when automatic registration of frame is resolved
        matrix.dml = [ self.ID,  ' = load(\" \", format=\"csv\")\n'] + matrix.dml
        if isinstance(self.eval_data, SUPPORTED_TYPES) and execute:
            matrix.script.input(self.ID, convertToMatrixBlock(matrix.sc, self.eval_data))
        elif execute:
            matrix.script.input(self.ID, self.toDF())
        return self

    def _register_as_output(self, execute):
        # TODO: Remove this when automatic registration of frame is resolved
        matrix.dml = matrix.dml + ['save(',  self.ID, ', \" \")\n']
        if execute:
            matrix.script.output(self.ID)

    def _visit(self, execute=True):
        """
        This function is called for two scenarios:
        1. For printing the PyDML script which has not yet been evaluated (execute=False). See '__repr__' method.
        2. Called as part of 'eval' method (execute=True). In this scenario, it builds the PyDML script by visiting itself
        and its child nodes. Also, it does appropriate registration as input or output that is required by MLContext.
        """
        if self.visited:
            return self
        self._mark_as_visited()
        if self.eval_data is not None:
            self._register_as_input(execute)
        elif self.op is not None:
            # Traverse the AST
            for m in self.op.inputs:
                m._visit(execute=execute)
            self.op._visit(execute=execute)
        else:
            raise Exception('Expected either op or data to be set')
        if self.eval_data is None and self.output:
            self._register_as_output(execute)
        return self

    def print_ast(self):
        """
        Please use m.print_ast() and/or  type `m` for debugging. Here is a sample session:

        >>> npm = np.ones((3,3))
        >>> m1 = sml.matrix(npm + 3)
        >>> m2 = sml.matrix(npm + 5)
        >>> m3 = m1 + m2
        >>> m3
        mVar2 = load(" ", format="csv")
        mVar1 = load(" ", format="csv")
        mVar3 = mVar1 + mVar2
        save(mVar3, " ")
        >>> m3.print_ast()
        - [mVar3] (op).
          - [mVar1] (data).
          - [mVar2] (data).
        """
        return self._print_ast(0)

    def _print_ast(self, numSpaces):
        head = ''.join([ ' ' ]*numSpaces + [ '- [', self.ID, '] ' ])
        if self.eval_data is not None:
            out = head + '(data).\n'
        elif self.op is not None:
            ret = [ head, '(op).\n' ]
            for m in self.op.inputs:
                ret = ret + [ m._print_ast(numSpaces + 2) ]
            out = ''.join(ret)
        else:
            raise ValueError('Either op or data needs to be set')
        if numSpaces == 0:
            print(out)
        else:
            return out

    def __repr__(self):
        """
        This function helps to debug matrix class and also examine the generated PyDML script
        """
        if self.eval_data is None:
            print('# This matrix (' + self.ID + ') is backed by below given PyDML script (which is not yet evaluated). To fetch the data of this matrix, invoke toNumPy() or toDF() or toPandas() methods.\n' + eval([self], execute=False))
        else:
            print('# This matrix (' + self.ID + ') is backed by ' + str(type(self.eval_data)) + '. To fetch the DataFrame or NumPy array, invoke toDF() or toNumPy() method respectively.')
        return ''

    ######################### NumPy related methods ######################################

    __array_priority__ = 10.2
    ndim = 2

    THROW_ARRAY_CONVERSION_ERROR = False

    def __array__(self, dtype=np.double):
        """
        As per NumPy from Python,
        This method is called to obtain an ndarray object when needed. You should always guarantee this returns an actual ndarray object.

        Using this method, you get back a ndarray object, and subsequent operations on the returned ndarray object will be singlenode.
        """
        if not isinstance(self.eval_data, SUPPORTED_TYPES):
            # Only warn if there is an unevaluated operation (which could potentially generate large matrix or if data is non-supported singlenode formats)
            import inspect
            frame,filename,line_number,function_name,lines,index = inspect.stack()[1]
            msg = 'Conversion from SystemML matrix to NumPy array (occurs in ' + str(filename) + ':' + str(line_number) + ' ' + function_name + ")"
            if matrix.THROW_ARRAY_CONVERSION_ERROR:
                raise Exception('[ERROR]:' + msg)
            else:
                print('[WARN]:' + msg)
        return np.array(self.toNumPy(), dtype)

    def astype(self, t):
        # TODO: Throw error if incorrect type
        return self

    def asfptype(self):
        return self

    def set_shape(self,shape):
        raise NotImplementedError('Reshaping is not implemented')

    def get_shape(self):
        if self._shape is None:
            lhsStr, inputs = _matricize(self, [])
            rlen_ID = _get_new_var_id()
            clen_ID = _get_new_var_id()
            multiline_dml = [rlen_ID, ' = ', lhsStr, '.shape(0)\n']
            multiline_dml = multiline_dml + [clen_ID, ' = ', lhsStr, '.shape(1)\n']
            multiline_dml = multiline_dml + [OUTPUT_ID, ' = full(0, rows=2, cols=1)\n']
            multiline_dml = multiline_dml + [ OUTPUT_ID, '[0,0] = ', rlen_ID, '\n' ]
            multiline_dml = multiline_dml + [ OUTPUT_ID, '[1,0] = ', clen_ID, '\n' ]
            ret = construct_intermediate_node(inputs, multiline_dml).toNumPy()
            self._shape = tuple(np.array(ret, dtype=int).flatten())
        return self._shape

    shape = property(fget=get_shape, fset=set_shape)

    def __numpy_ufunc__(self, func, method, pos, inputs, **kwargs):
        """
        This function enables systemml matrix to be compatible with NumPy's ufuncs.

        Parameters
        ----------
        func:  ufunc object that was called.
        method: string indicating which Ufunc method was called (one of "__call__", "reduce", "reduceat", "accumulate", "outer", "inner").
        pos: index of self in inputs.
        inputs:  tuple of the input arguments to the ufunc
        kwargs: dictionary containing the optional input arguments of the ufunc.
        """
        if method != '__call__' or kwargs:
            return NotImplemented
        if func in matrix._numpy_to_systeml_mapping:
            fn = matrix._numpy_to_systeml_mapping[func]
        else:
            return NotImplemented
        if len(inputs) == 2:
            return fn(inputs[0], inputs[1])
        elif  len(inputs) == 1:
            return fn(inputs[0])
        else:
            raise ValueError('Unsupported number of inputs')

    def hstack(self, other):
        """
        Stack matrices horizontally (column wise). Invokes cbind internally.
        """
        return binaryMatrixFunction(self, other, 'cbind')

    def vstack(self, other):
        """
        Stack matrices vertically (row wise). Invokes rbind internally.
        """
        return binaryMatrixFunction(self, other, 'rbind')

    ######################### Arithmetic operators ######################################

    def negative(self):
        lhsStr, inputs = _matricize(self, [])
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = -', lhsStr, '\n'])

    def remainder(self, other):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        rhsStr, inputs = _matricize(other, inputs)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = floor(', lhsStr, '/', rhsStr, ') * ', rhsStr, '\n'])

    def ldexp(self, other):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        rhsStr, inputs = _matricize(other, inputs)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = ', lhsStr, '* (2**', rhsStr, ')\n'])

    def mod(self, other):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        rhsStr, inputs = _matricize(other, inputs)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = ', lhsStr, ' - floor(', lhsStr, '/', rhsStr, ') * ', rhsStr, '\n'])

    def logaddexp(self, other):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        rhsStr, inputs = _matricize(other, inputs)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = log(exp(', lhsStr, ') + exp(', rhsStr, '))\n'])

    def logaddexp2(self, other):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        rhsStr, inputs = _matricize(other, inputs)
        opStr =  _log_base('2**' + lhsStr + '2**' + rhsStr, 2)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = ', opStr, '\n'])

    def log1p(self):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = log(1 + ', lhsStr, ')\n'])

    def exp(self):
        return unaryMatrixFunction(self, 'exp')

    def exp2(self):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = 2**', lhsStr, '\n'])

    def square(self):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = ', lhsStr, '**2\n'])

    def reciprocal(self):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = 1/', lhsStr, '\n'])

    def expm1(self):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = exp(', lhsStr, ') - 1\n'])

    def ones_like(self):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        rlen = lhsStr + '.shape(axis=0)'
        clen = lhsStr + '.shape(axis=1)'
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = full(1, rows=', rlen, ', cols=', clen, ')\n'])

    def zeros_like(self):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        rlen = lhsStr + '.shape(axis=0)'
        clen = lhsStr + '.shape(axis=1)'
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = full(0, rows=', rlen, ', cols=', clen, ')\n'])

    def log2(self):
        return self.log(2)

    def log10(self):
        return self.log(10)

    def log(self, y=None):
        if y is None:
            return unaryMatrixFunction(self, 'log')
        else:
            return binaryMatrixFunction(self, y, 'log')

    def abs(self):
        return unaryMatrixFunction(self, 'abs')

    def sqrt(self):
        return unaryMatrixFunction(self, 'sqrt')

    def round(self):
        return unaryMatrixFunction(self, 'round')

    def floor(self):
        return unaryMatrixFunction(self, 'floor')

    def ceil(self):
        return unaryMatrixFunction(self, 'ceil')

    def sin(self):
        return unaryMatrixFunction(self, 'sin')

    def cos(self):
        return unaryMatrixFunction(self, 'cos')

    def tan(self):
        return unaryMatrixFunction(self, 'tan')

    def arcsin(self):
        return self.asin()

    def arccos(self):
        return self.acos()

    def arctan(self):
        return self.atan()

    def asin(self):
        return unaryMatrixFunction(self, 'asin')

    def acos(self):
        return unaryMatrixFunction(self, 'acos')

    def atan(self):
        return unaryMatrixFunction(self, 'atan')

    def rad2deg(self):
        """
        Convert angles from radians to degrees.
        """
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        # 180/pi = 57.2957795131
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = ', lhsStr, '*57.2957795131\n'])

    def deg2rad(self):
        """
        Convert angles from degrees to radians.
        """
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        # pi/180 = 0.01745329251
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = ', lhsStr, '*0.01745329251\n'])

    def sign(self):
        return unaryMatrixFunction(self, 'sign')

    def __add__(self, other):
        return binary_op(self, other, ' + ')

    def __sub__(self, other):
        return binary_op(self, other, ' - ')

    def __mul__(self, other):
        return binary_op(self, other, ' * ')

    def __floordiv__(self, other):
        return binary_op(self, other, ' // ')

    def __div__(self, other):
        """
        Performs division (Python 2 way).
        """
        return binary_op(self, other, ' / ')

    def __truediv__(self, other):
        """
        Performs division (Python 3 way).
        """
        return binary_op(self, other, ' / ')

    def __mod__(self, other):
        return binary_op(self, other, ' % ')

    def __pow__(self, other):
        return binary_op(self, other, ' ** ')

    def __radd__(self, other):
        return binary_op(other, self, ' + ')

    def __rsub__(self, other):
        return binary_op(other, self, ' - ')

    def __rmul__(self, other):
        return binary_op(other, self, ' * ')

    def __rfloordiv__(self, other):
        return binary_op(other, self, ' // ')

    def __rdiv__(self, other):
        return binary_op(other, self, ' / ')

    def __rmod__(self, other):
        return binary_op(other, self, ' % ')

    def __rpow__(self, other):
        return binary_op(other, self, ' ** ')

    def dot(self, other):
        """
        Numpy way of performing matrix multiplication
        """
        return binaryMatrixFunction(self, other, 'dot')

    def __matmul__(self, other):
        """
        Performs matrix multiplication (infix operator: @). See PEP 465)
        """
        return binaryMatrixFunction(self, other, 'dot')

    ######################### Relational/Boolean operators ######################################

    def __lt__(self, other):
        return binary_op(self, other, ' < ')

    def __le__(self, other):
        return binary_op(self, other, ' <= ')

    def __gt__(self, other):
        return binary_op(self, other, ' > ')

    def __ge__(self, other):
        return binary_op(self, other,' >= ')

    def __eq__(self, other):
        return binary_op(self, other, ' == ')

    def __ne__(self, other):
        return binary_op(self, other, ' != ')

    # TODO: Cast the output back into scalar and return boolean results
    def __and__(self, other):
        return binary_op(other, self, ' & ')

    def __or__(self, other):
        return binary_op(other, self, ' | ')

    def logical_not(self):
        inputs = []
        lhsStr, inputs = _matricize(self, inputs)
        return construct_intermediate_node(inputs, [OUTPUT_ID, ' = !', lhsStr, '\n'])

    def remove_empty(self, axis=None):
        """
        Removes all empty rows or columns from the input matrix target X according to specified axis.

        Parameters
        ----------
        axis : int (0 or 1)
        """
        if axis is None:
            raise ValueError('axis is a mandatory argument for remove_empty')
        if axis == 0:
            return self._parameterized_helper_fn(self, 'removeEmpty',  { 'target':self, 'margin':'rows' })
        elif axis == 1:
            return self._parameterized_helper_fn(self, 'removeEmpty',  { 'target':self, 'margin':'cols' })
        else:
            raise ValueError('axis for remove_empty needs to be either 0 or 1.')

    def replace(self, pattern=None, replacement=None):
        """
        Removes all empty rows or columns from the input matrix target X according to specified axis.

        Parameters
        ----------
        pattern : float or int
        replacement : float or int
        """
        if pattern is None or not isinstance(pattern, (float, int)):
            raise ValueError('pattern should be of type float or int')
        if replacement is None or not isinstance(replacement, (float, int)):
            raise ValueError('replacement should be of type float or int')
        return self._parameterized_helper_fn(self, 'replace',  { 'target':self, 'pattern':pattern, 'replacement':replacement })

    def _parameterized_helper_fn(self, fnName, **kwargs):
        """
        Helper to invoke parameterized builtin function
        """
        dml_script = ''
        lhsStr, inputs = _matricize(self, [])
        dml_script = [OUTPUT_ID, ' = ', fnName, '(', lhsStr ]
        first_arg = True
        for key in kwargs:
            if first_arg:
                first_arg = False
            else:
                dml_script = dml_script + [ ', ' ]
            v = kwargs[key]
            if isinstance(v, str):
                dml_script = dml_script + [key, '=\"', v, '\"' ]
            elif isinstance(v, matrix):
                dml_script = dml_script + [key, '=', v.ID]
            else:
                dml_script = dml_script + [key, '=', str(v) ]
        dml_script = dml_script + [ ')\n' ]
        return construct_intermediate_node(inputs, dml_script)

    ######################### Aggregation functions ######################################

    def prod(self):
        """
        Return the product of all cells in matrix
        """
        return self._aggFn('prod', None)

    def sum(self, axis=None):
        """
        Compute the sum along the specified axis

        Parameters
        ----------
        axis : int, optional
        """
        return self._aggFn('sum', axis)

    def mean(self, axis=None):
        """
        Compute the arithmetic mean along the specified axis

        Parameters
        ----------
        axis : int, optional
        """
        return self._aggFn('mean', axis)

    def var(self, axis=None):
        """
        Compute the variance along the specified axis.
        We assume that delta degree of freedom is 1 (unlike NumPy which assumes ddof=0).

        Parameters
        ----------
        axis : int, optional
        """
        return self._aggFn('var', axis)

    def moment(self, moment=1, axis=None):
        """
        Calculates the nth moment about the mean

        Parameters
        ----------
        moment : int
            can be 1, 2, 3 or 4
        axis : int, optional
        """
        if moment == 1:
            return self.mean(axis)
        elif moment == 2:
            return self.var(axis)
        elif moment == 3 or moment == 4:
            return self._moment_helper(moment, axis)
        else:
            raise ValueError('The specified moment is not supported:' + str(moment))

    def _moment_helper(self, k, axis=0):
        dml_script = ''
        lhsStr, inputs = _matricize(self, [])
        dml_script = [OUTPUT_ID, ' = moment(', lhsStr, ', ', str(k), ')\n' ]
        dml_script = [OUTPUT_ID, ' = moment(', lhsStr, ', ', str(k), ')\n' ]
        if axis is None:
            dml_script = [OUTPUT_ID, ' = moment(full(', lhsStr, ', rows=length(', lhsStr, '), cols=1), ', str(k), ')\n' ]
        elif axis == 0:
            dml_script = [OUTPUT_ID, ' = full(0, rows=nrow(', lhsStr, '), cols=1)\n' ]
            dml_script = dml_script + [ 'parfor(i in 1:nrow(', lhsStr, '), check=0):\n' ]
            dml_script = dml_script + [ '\t', OUTPUT_ID, '[i-1, 0] = moment(full(', lhsStr, '[i-1,], rows=ncol(', lhsStr, '), cols=1), ', str(k), ')\n\n' ]
        elif axis == 1:
            dml_script = [OUTPUT_ID, ' = full(0, rows=1, cols=ncol(', lhsStr, '))\n' ]
            dml_script = dml_script + [ 'parfor(i in 1:ncol(', lhsStr, '), check=0):\n' ]
            dml_script = dml_script + [ '\t', OUTPUT_ID, '[0, i-1] = moment(', lhsStr, '[,i-1], ', str(k), ')\n\n' ]
        else:
            raise ValueError('Incorrect axis:' + axis)
        return construct_intermediate_node(inputs, dml_script)

    def sd(self, axis=None):
        """
        Compute the standard deviation along the specified axis

        Parameters
        ----------
        axis : int, optional
        """
        return self._aggFn('sd', axis)

    def max(self, other=None, axis=None):
        """
        Compute the maximum value along the specified axis

        Parameters
        ----------
        other: matrix or numpy array (& other supported types) or scalar
        axis : int, optional
        """
        if other is not None and axis is not None:
            raise ValueError('Both axis and other cannot be not None')
        elif other is None and axis is not None:
            return self._aggFn('max', axis)
        else:
            return binaryMatrixFunction(self, other, 'max')

    def min(self, other=None, axis=None):
        """
        Compute the minimum value along the specified axis

        Parameters
        ----------
        other: matrix or numpy array (& other supported types) or scalar
        axis : int, optional
        """
        if other is not None and axis is not None:
            raise ValueError('Both axis and other cannot be not None')
        elif other is None and axis is not None:
            return self._aggFn('min', axis)
        else:
            return binaryMatrixFunction(self, other, 'min')

    def argmin(self, axis=None):
        """
        Returns the indices of the minimum values along an axis.

        Parameters
        ----------
        axis : int, optional  (only axis=1, i.e. rowIndexMax is supported in this version)
        """
        return self._aggFn('argmin', axis)

    def argmax(self, axis=None):
        """
        Returns the indices of the maximum values along an axis.

        Parameters
        ----------
        axis : int, optional (only axis=1, i.e. rowIndexMax is supported in this version)
        """
        return self._aggFn('argmax', axis)

    def cumsum(self, axis=None):
        """
        Returns the indices of the maximum values along an axis.

        Parameters
        ----------
        axis : int, optional (only axis=0, i.e. cumsum along the rows is supported in this version)
        """
        return self._aggFn('cumsum', axis)

    def transpose(self):
        """
        Transposes the matrix.
        """
        return self._aggFn('transpose', None)

    def trace(self):
        """
        Return the sum of the cells of the main diagonal square matrix
        """
        return self._aggFn('trace', None)

    def _aggFn(self, fnName, axis):
        """
        Common function that is called for functions that have axis as parameter.
        """
        dml_script = ''
        lhsStr, inputs = _matricize(self, [])
        if axis is None:
            dml_script = [OUTPUT_ID, ' = ', fnName, '(', lhsStr, ')\n']
        else:
            dml_script = [OUTPUT_ID, ' = ', fnName, '(', lhsStr, ', axis=', str(axis) ,')\n']
        return construct_intermediate_node(inputs, dml_script)

    ######################### Indexing operators ######################################

    def __getitem__(self, index):
        """
        Implements evaluation of right indexing operations such as m[1,1], m[0:1,], m[:, 0:1]
        """
        return construct_intermediate_node([self], [OUTPUT_ID, ' = ', self.ID ] + getIndexingDML(index) + [ '\n' ])

    # Performs deep copy if the matrix is backed by data
    def _prepareForInPlaceUpdate(self):
        temp = matrix(self.eval_data, op=self.op)
        for op in self.referenced:
            op.inputs = [temp if x.ID==self.ID else x for x in op.inputs]
        self.ID, temp.ID = temp.ID, self.ID # Copy even the IDs as the IDs might be used to create DML
        self.op = DMLOp([temp], dml=[self.ID, " = ", temp.ID])
        self.eval_data = None
        temp.referenced = self.referenced + [ self.op ]
        self.referenced = []

    def __setitem__(self, index, value):
        """
        Implements evaluation of left indexing operations such as m[1,1]=2
        """
        self._prepareForInPlaceUpdate()
        if isinstance(value, matrix) or isinstance(value, DMLOp):
            self.op.inputs = self.op.inputs + [ value ]
        if isinstance(value, matrix):
            value.referenced = value.referenced + [ self.op ]
        self.op.dml = self.op.dml + [ '\n', self.ID ] + getIndexingDML(index) + [ ' = ',  getValue(value), '\n']

    # Not implemented: conj, hyperbolic/inverse-hyperbolic functions(i.e. sinh, arcsinh, cosh, ...), bitwise operator, xor operator, isreal, iscomplex, isfinite, isinf, isnan, copysign, nextafter, modf, frexp, trunc
    _numpy_to_systeml_mapping = {np.add: __add__, np.subtract: __sub__, np.multiply: __mul__, np.divide: __div__, np.logaddexp: logaddexp, np.true_divide: __truediv__, np.floor_divide: __floordiv__, np.negative: negative, np.power: __pow__, np.remainder: remainder, np.mod: mod, np.fmod: __mod__, np.absolute: abs, np.rint: round, np.sign: sign, np.exp: exp, np.exp2: exp2, np.log: log, np.log2: log2, np.log10: log10, np.expm1: expm1, np.log1p: log1p, np.sqrt: sqrt, np.square: square, np.reciprocal: reciprocal, np.ones_like: ones_like, np.zeros_like: zeros_like, np.sin: sin, np.cos: cos, np.tan: tan, np.arcsin: arcsin, np.arccos: arccos, np.arctan: arctan, np.deg2rad: deg2rad, np.rad2deg: rad2deg, np.greater: __gt__, np.greater_equal: __ge__, np.less: __lt__, np.less_equal: __le__, np.not_equal: __ne__, np.equal: __eq__, np.logical_not: logical_not, np.logical_and: __and__, np.logical_or: __or__, np.maximum: max, np.minimum: min, np.signbit: sign, np.ldexp: ldexp, np.dot:dot}

ext: py

max_stars_repo_name: ckadner/incubator-systemml

max_issues_repo_name: ckadner/incubator-systemml

max_forks_repo_name: ckadner/incubator-systemml

reason: import numpy,from scipy


text:
# Copyright 2019 Jij Inc.
# 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

import numpy as np
import openjij
from openjij.sampler import BaseSampler
from openjij.utils.decorator import deprecated_alias
from openjij.utils.graph_utils import qubo_to_ising
from openjij.model import BinaryHigherOrderModel
from .hubo_simulated_annealing import hubo_sa_sampling
from .hubo_simulated_annealing import default_schedule
import cxxjij

"""
This module contains Simulated Annealing sampler.
"""

class SASampler(BaseSampler):
    """Sampler with Simulated Annealing (SA).

    Args:
        beta_min (float):
            Minmum beta (inverse temperature).
            You can overwrite in methods .sample_*.

        beta_max (float):
            Maximum beta (inverse temperature).
            You can overwrite in methods .sample_*.

        num_reads (int):
            number of sampling (algorithm) runs. defaults None.
            You can overwrite in methods .sample_*.

        num_sweeps (int):
            number of MonteCarlo steps during SA. defaults None.
            You can overwrite in methods .sample_*.

        schedule_info (dict):
            Information about an annealing schedule.

    Raises:
        ValueError: If schedules or variables violate as below.
        - not list or numpy.array.
        - not list of tuple (beta : float, step_length : int).
        - beta is less than zero.

    """

    @property
    def parameters(self):
        return {
            'beta_min': ['parameters'],
            'beta_max': ['parameters'],
        }

    def __init__(self,
                 beta_min=None, beta_max=None,
                 num_sweeps=1000, schedule=None,
                 num_reads=1):

        self.beta_min = beta_min
        self.beta_max = beta_max
        self.num_sweeps = num_sweeps
        self.schedule = schedule
        self.num_reads = num_reads
        self._schedule_setting = {
            'beta_min': beta_min,
            'beta_max': beta_max,
            'num_sweeps': num_sweeps,
            'num_reads': num_reads,
        }

        self._make_system = {
            'singlespinflip': cxxjij.system.make_classical_ising,
            'swendsenwang': cxxjij.system.make_classical_ising
        }
        self._algorithm = {
            'singlespinflip': cxxjij.algorithm.Algorithm_SingleSpinFlip_run,
            'swendsenwang': cxxjij.algorithm.Algorithm_SwendsenWang_run
        }

    def _convert_validation_schedule(self, schedule):
        """Checks if the schedule is valid and returns cxxjij schedule
        """
        if not isinstance(schedule, (list, np.array)):
            raise ValueError("schedule should be list or numpy.array")

        if isinstance(schedule[0], cxxjij.utility.ClassicalSchedule):
            return schedule

        if len(schedule[0]) != 2:
            raise ValueError(
                "schedule is list of tuple or list (beta : float, step_length : int)")

        # schedule validation  0 <= beta
        beta = np.array(schedule).T[0]
        if not np.all(0 <= beta):
            raise ValueError("schedule beta range is '0 <= beta'.")

        # convert to cxxjij.utility.ClassicalSchedule
        cxxjij_schedule = []
        for beta, step_length in schedule:
            _schedule = cxxjij.utility.ClassicalSchedule()
            _schedule.one_mc_step = step_length
            _schedule.updater_parameter.beta = beta
            cxxjij_schedule.append(_schedule)

        return cxxjij_schedule

    def sample_ising(self, h, J, beta_min=None, beta_max=None,
                     num_sweeps=None, num_reads=1, schedule=None,
                     initial_state=None, updater='single spin flip',
                     sparse=False,
                     reinitialize_state=True, seed=None,
                     ):
        """sample Ising model.

        Args:
            h (dict): linear biases
            J (dict): quadratic biases
            beta_min (float): minimal value of inverse temperature
            beta_max (float): maximum value of inverse temperature
            num_sweeps (int): number of sweeps
            num_reads (int): number of reads
            schedule (list): list of inverse temperature
            initial_state (dict): initial state
            updater(str): updater algorithm
            reinitialize_state (bool): if true reinitialize state for each run
            seed (int): seed for Monte Carlo algorithm
        Returns:
            :class:`openjij.sampler.response.Response`: results

        Examples:

            for Ising case::

                >>> h = {0: -1, 1: -1, 2: 1, 3: 1}
                >>> J = {(0, 1): -1, (3, 4): -1}
                >>> sampler = oj.SASampler()
                >>> res = sampler.sample_ising(h, J)

            for QUBO case::

                >>> Q = {(0, 0): -1, (1, 1): -1, (2, 2): 1, (3, 3): 1, (4, 4): 1, (0, 1): -1, (3, 4): 1}
                >>> sampler = oj.SASampler()
                >>> res = sampler.sample_qubo(Q)

        """

        model = openjij.BinaryQuadraticModel(
            linear=h, quadratic=J, var_type='SPIN'
        )
        return self._sampling(model, beta_min, beta_max,
                              num_sweeps, num_reads, schedule,
                              initial_state, updater,
                              sparse=sparse,
                              reinitialize_state=reinitialize_state,
                              seed=seed)

    def _sampling(self, model, beta_min=None, beta_max=None,
                     num_sweeps=None, num_reads=1, schedule=None,
                     initial_state=None, updater='single spin flip',
                     sparse=False,
                     reinitialize_state=True, seed=None, structure=None,
                     ):
        """sampling by using specified model
        Args:
            model (openjij.BinaryQuadraticModel): BinaryQuadraticModel
            beta_min (float): minimal value of inverse temperature
            beta_max (float): maximum value of inverse temperature
            num_sweeps (int): number of sweeps
            num_reads (int): number of reads
            schedule (list): list of inverse temperature
            initial_state (dict): initial state
            updater(str): updater algorithm
            reinitialize_state (bool): if true reinitialize state for each run
            seed (int): seed for Monte Carlo algorithm
            structure (dict): specify the structure.
            This argument is necessary if the model has a specific structure (e.g. Chimera graph) and the updater algorithm is structure-dependent.
            structure must have two types of keys, namely "size" which shows the total size of spins and "dict" which is the map from model index (elements in model.indices) to the number.
        Returns:
            :class:`openjij.sampler.response.Response`: results
        """
        _updater_name = updater.lower().replace('_', '').replace(' ', '')
        # swendsen wang algorithm runs only on sparse ising graphs.
        if _updater_name == 'swendsenwang' or sparse:
            ising_graph = model.get_cxxjij_ising_graph(sparse=True)
        else:
            ising_graph = model.get_cxxjij_ising_graph()

        self._setting_overwrite(
            beta_min=beta_min, beta_max=beta_max,
            num_sweeps=num_sweeps, num_reads=num_reads
        )

        # set annealing schedule -------------------------------
        if schedule or self.schedule:
            self._schedule = self._convert_validation_schedule(
                schedule if schedule else self.schedule
            )
            self.schedule_info = {'schedule': 'custom schedule'}
        else:
            self._schedule, beta_range = geometric_ising_beta_schedule(
                model=model,
                beta_max=self._schedule_setting['beta_max'],
                beta_min=self._schedule_setting['beta_min'],
                num_sweeps=self._schedule_setting['num_sweeps']
            )
            self.schedule_info = {
                'beta_max': beta_range[0],
                'beta_min': beta_range[1],
                'num_sweeps': self._schedule_setting['num_sweeps']
            }
        # ------------------------------- set annealing schedule

        # make init state generator --------------------------------
        if initial_state is None:
            def _generate_init_state(): return ising_graph.gen_spin(seed) if seed != None else ising_graph.gen_spin()
        else:
            if isinstance(initial_state, dict):
                initial_state = [initial_state[k] for k in model.indices]
            _init_state = np.array(initial_state)

            if structure == None:
                # validate initial_state size
                if len(initial_state) != ising_graph.size():
                    raise ValueError(
                        "the size of the initial state should be {}"
                        .format(ising_graph.size()))
            else:
                # resize _initial_state
                temp_state = [1]*int(structure['size'])
                for k,ind in enumerate(model.indices):
                    temp_state[structure['dict'][ind]] = _init_state[k]
                _init_state = temp_state

            def _generate_init_state(): return np.array(_init_state)
        # -------------------------------- make init state generator

        # choose updater -------------------------------------------
        _updater_name = updater.lower().replace('_', '').replace(' ', '')
        if _updater_name not in self._make_system:
            raise ValueError('updater is one of "single spin flip or swendsen wang"')
        algorithm = self._algorithm[_updater_name]
        sa_system = self._make_system[_updater_name](_generate_init_state(), ising_graph)
        # ------------------------------------------- choose updater
        response = self._cxxjij_sampling(
            model, _generate_init_state,
            algorithm, sa_system,
            reinitialize_state, seed, structure
        )

        response.info['schedule'] = self.schedule_info

        return response

    def sample_hubo(self, interactions: list, var_type,
                    beta_min=None, beta_max=None, schedule=None,
                    num_sweeps=100, num_reads=1,
                    init_state=None, seed=None):
        """sampling from higher order unconstrainted binary optimization.

        Args:
            interactions (list): ordered by degree of interaction.
            var_type (str, openjij.VarType): "SPIN" or "BINARY"
            beta_min (float, optional): Minimum beta (initial inverse temperature). Defaults to None.
            beta_max (float, optional): Maximum beta (final inverse temperature). Defaults to None.
            schedule (list, optional): schedule list. Defaults to None.
            num_sweeps (int, optional): number of sweeps. Defaults to 100.
            num_reads (int, optional): number of reads. Defaults to 1.
            init_state (list, optional): initial state. Defaults to None.
            seed (int, optional): [description]. Defaults to None.

        Returns:
            :class:`openjij.sampler.response.Response`: results

        Examples::

            >>> sampler = oj.SASampler()
            >>> h = {0: -1}
            >>> J = {(0, 1): -1}
            >>> K = {(0, 1, 2): 1}
            >>> response = sampler.sample_hubo([h, J, K], var_type="SPIN")

        """

        self._setting_overwrite(
            beta_min=beta_min, beta_max=beta_max,
            num_sweeps=num_sweeps, num_reads=num_reads
        )

        bhom = BinaryHigherOrderModel(interactions)

        if schedule or self.schedule:
            self._schedule = self._convert_validation_schedule(
                schedule if schedule else self.schedule
            )
            self.schedule_info = {'schedule': 'custom schedule'}
        else:
            schedule = default_schedule(
                bhom,
                beta_min=self._schedule_setting['beta_min'],
                beta_max=self._schedule_setting['beta_max'],
                num_sweeps=self._schedule_setting['num_sweeps'])
            self.schedule_info = {
                'beta_max': schedule[-1][0],
                'beta_min': schedule[0][0],
                'num_sweeps': self._schedule_setting['num_sweeps']
            }

        response = hubo_sa_sampling(
            bhom, var_type,
            schedule=schedule, schedule_info=self.schedule_info,
            num_sweeps=self._schedule_setting['num_sweeps'],
            num_reads=self._schedule_setting['num_reads'],
            init_state=init_state, seed=seed
        )
        return response

def geometric_ising_beta_schedule(model: openjij.model.BinaryQuadraticModel,
                                  beta_max=None, beta_min=None,
                                  num_sweeps=1000):
    """make geometric cooling beta schedule

    Args:
        model (openjij.BinaryQuadraticModel)
        beta_max (float, optional): [description]. Defaults to None.
        beta_min (float, optional): [description]. Defaults to None.
        num_sweeps (int, optional): [description]. Defaults to 1000.
    Returns:
        list of cxxjij.utility.ClassicalSchedule, list of beta range [max, min]
    """
    if beta_min is None or beta_max is None:
        # generate Ising matrix
        ising_interaction = model.interaction_matrix()
        if (model.vartype == openjij.BINARY):
            # convert to ising matrix
            qubo_to_ising(ising_interaction)

        ising_interaction = np.abs(ising_interaction)

        #automatical setting of min, max delta energy
        abs_bias = np.sum(ising_interaction, axis=1)

        min_delta_energy = np.min(ising_interaction[ising_interaction > 0])
        max_delta_energy = np.max(abs_bias[abs_bias > 0])

    # TODO: More optimal schedule ?

    beta_min = np.log(2) / max_delta_energy if beta_min is None else beta_min
    beta_max = np.log(100) / min_delta_energy if beta_max is None else beta_max

    num_sweeps_per_beta = max(1, num_sweeps // 1000)

    # set schedule to cxxjij
    schedule = cxxjij.utility.make_classical_schedule_list(
        beta_min=beta_min, beta_max=beta_max,
        one_mc_step=num_sweeps_per_beta,
        num_call_updater=num_sweeps//num_sweeps_per_beta
    )

    return schedule, [beta_max, beta_min]

reason: import numpy


text:
__author__ = 'sibirrer'

import numpy as np
import numpy.testing as npt
import pytest
from lenstronomy.LensModel.Profiles.curved_arc_sis_mst import CurvedArcSISMST
from lenstronomy.LensModel.Profiles.curved_arc_const import CurvedArcConstMST, CurvedArcConst
from lenstronomy.Util import util

class TestCurvedArcConst(object):
    """
    tests the source model routines
    """
    def setup(self):
        self.arc_sis = CurvedArcSISMST()
        self.arc_const = CurvedArcConst()

    def test_function(self):
        kwargs_arc = {'tangential_stretch': 5,
                      #'radial_stretch': 1,
                      'curvature': 1. / 10,
                      'direction': 0,
                      'center_x': 0,
                      'center_y': 0
                      }
        npt.assert_raises(Exception, self.arc_const.function, 0., 0., **kwargs_arc)

    def test_derivatives(self):
        kwargs_arc = {'tangential_stretch': 3,
                      #'radial_stretch': 1.,
                      'curvature': 0.8,
                      'direction': 0,
                      'center_x': 0,
                      'center_y': 0
                      }

        kwargs_arc_sis = {'tangential_stretch': 3,
                      'radial_stretch': 1.,
                      'curvature': 0.8,
                      'direction': 0,
                      'center_x': 0,
                      'center_y': 0
                      }
        x, y = util.make_grid(numPix=100, deltapix=0.01)
        f_x_sis, f_y_sis = self.arc_sis.derivatives(x, y, **kwargs_arc_sis)
        beta_x_sis = x - f_x_sis
        beta_y_sis = y - f_y_sis
        f_x_const, f_y_const = self.arc_const.derivatives(x, y, **kwargs_arc)
        beta_x_const = x - f_x_const
        beta_y_const = y - f_y_const

        from lenstronomy.LightModel.light_model import LightModel
        gauss = LightModel(['GAUSSIAN'])
        kwargs_source = [{'amp': 1, 'sigma': 0.05, 'center_x': 0, 'center_y': 0}]
        flux_sis = gauss.surface_brightness(beta_x_sis, beta_y_sis, kwargs_source)
        flux_const = gauss.surface_brightness(beta_x_const, beta_y_const, kwargs_source)

        npt.assert_almost_equal((flux_const - flux_sis) / np.max(flux_const), 0, decimal=2)

        # check for stability outside the defined bounds of curvature
        f_x_const, f_y_const = self.arc_const.derivatives(x=0, y=1000, **kwargs_arc)
        npt.assert_almost_equal(f_x_const, 0)
        npt.assert_almost_equal(f_y_const, 0)

    def test_hessian(self):
        kwargs_arc = {'tangential_stretch': 5,
                      #'radial_stretch': 1,
                      'curvature': 1. / 10,
                      'direction': 0.5,
                      'center_x': 0,
                      'center_y': 0
                      }
        x, y = 0., 1.
        f_xx, f_xy, f_yx, f_yy = self.arc_const.hessian(x, y, **kwargs_arc)

        alpha_ra, alpha_dec = self.arc_const.derivatives(x, y, **kwargs_arc)
        diff = 0.0000001
        alpha_ra_dx, alpha_dec_dx = self.arc_const.derivatives(x + diff, y, **kwargs_arc)
        alpha_ra_dy, alpha_dec_dy = self.arc_const.derivatives(x, y + diff, **kwargs_arc)

        f_xx_num = (alpha_ra_dx - alpha_ra) / diff
        f_xy_num = (alpha_ra_dy - alpha_ra) / diff
        f_yx_num = (alpha_dec_dx - alpha_dec) / diff
        f_yy_num = (alpha_dec_dy - alpha_dec) / diff
        print(f_xx, f_xx_num)
        print(f_xy, f_xy_num)
        print(f_yx, f_yx_num)
        print(f_yy, f_yy_num)

        npt.assert_almost_equal(f_xx, f_xx_num)
        npt.assert_almost_equal(f_xy, f_xy_num)
        npt.assert_almost_equal(f_yx, f_yx_num)
        npt.assert_almost_equal(f_yy, f_yy_num)

class TestCurvedArcConstMST(object):

    def setup(self):
        self.arc_const = CurvedArcConstMST()

    def test_function(self):
        kwargs_arc = {'tangential_stretch': 5,
                      'radial_stretch': 1,
                      'curvature': 1. / 10,
                      'direction': 0,
                      'center_x': 0,
                      'center_y': 0
                      }
        npt.assert_raises(Exception, self.arc_const.function, 0., 0., **kwargs_arc)

    def test_hessian(self):
        kwargs_arc = {'tangential_stretch': 5,
                      'radial_stretch': 1,
                      'curvature': 1. / 10,
                      'direction': 0,
                      'center_x': 0,
                      'center_y': 0
                      }
        #npt.assert_raises(Exception, self.arc_const.hessian, 0., 0., **kwargs_arc)

if __name__ == '__main__':
    pytest.main("-k TestLensModel")


text:
# -*- coding: utf-8 -*-
from autograd import elementwise_grad as egrad  # for functions that vectorize over inputs
import matplotlib.pyplot as plt
import autograd.numpy as np  # Thinly-wrapped numpy
from autograd import grad    # The only autograd function you may ever need
import numpy
import sympy as sp
#最速下降法

THRESHOLD = 0.01
def func(x):                 # Define a function
    return  3*x * x + 16 * x - 9

grad_func = grad(func)       # Obtain its gradient function
x = -3.0
count = 0
while True:
    grad_val = grad_func(x)
    print 'outcome is  ', func(x)
    print x, ' **x'
    print grad_val, type(grad_val)
    if numpy.abs(grad_val) > THRESHOLD:
        d = -grad_val
        b = sp.Symbol("b")
        _x = sp.Symbol("_x")
        #得到一阶导数形式
        form =  sp.diff(func(_x), _x)
        form = form.replace('_x', '%f - b * %f' % (x, grad_val))
        #form = sp.sympify(form)
        #gamma 值,sp.S()转换字符串为函数形式，用于求解
        gamma = sp.solve(form, b)[0]
        #print 'gamma ', gamma
        x = x + gamma * d
        x = float(x)
        count = count + 1
        print '----------the ', count, 'th iteration-------------'
    else:
        print  ' Final location is ', final_X_Y, ' and the minimum is ', z.subs({a:final_X_Y[0], b:final_X_Y[1]})
        break

#linspacke 用于生成等差数列，第三个参数表示数列中数的个数
x = np.linspace(-7, 7, 200)
plt.plot(x, func(x))
plt.plot(x, egrad(func)(x))
'''x, egrad(tanh)(x),                                     # first  derivative
x, egrad(egrad(tanh))(x),                              # second derivative
x, egrad(egrad(egrad(tanh)))(x),                       # third  derivative
x, egrad(egrad(egrad(egrad(tanh))))(x),                # fourth derivative
x, egrad(egrad(egrad(egrad(egrad(tanh)))))(x),         # fifth  derivative
x, egrad(egrad(egrad(egrad(egrad(egrad(tanh))))))(x))'''
plt.show()

max_stars_repo_name: chuchienshu/Optimization-Assignments

max_issues_repo_name: chuchienshu/Optimization-Assignments

max_forks_repo_name: chuchienshu/Optimization-Assignments

reason: import numpy,import sympy


text:
# -*- coding: utf-8 -*-
"""
@author: Raluca Sandu
"""
from mpl_toolkits.mplot3d import axes3d
import cvxpy as cp
import matplotlib.pyplot as plt
import numpy as np
# import sklearn.datasets
from scipy.spatial import ConvexHull

def random_point_ellipsoid(a, b, c, x0, y0, z0):
    """Generate a random point on an ellipsoid defined by a,b,c"""
    u = np.random.rand()
    v = np.random.rand()
    theta = u * 2.0 * np.pi
    phi = np.arccos(2.0 * v - 1.0)
    sinTheta = np.sin(theta);
    cosTheta = np.cos(theta);
    sinPhi = np.sin(phi);
    cosPhi = np.cos(phi);
    rx = a * sinPhi * cosTheta;
    ry = b * sinPhi * sinTheta;
    rz = c * cosPhi;
    return rx, ry, rz

def random_point_ellipse(W, d):
    # random angle
    alpha = 2 * np.pi * np.random.random()
    # vector on that angle
    pt = np.array([np.cos(alpha), np.sin(alpha)])
    # Ellipsoidize it
    return W @ pt + d

def GetRandom(dims, Npts):
    if dims == 2:
        W = sklearn.datasets.make_spd_matrix(2)
        d = np.array([2, 3])
        points = np.array([random_point_ellipse(W, d) for i in range(Npts)])
    elif dims == 3:
        points = np.array([random_point_ellipsoid(3, 5, 7, 2, 3, 3) for i in range(Npts)])
    else:
        raise Exception("dims must be 2 or 3!")
    noise = np.random.multivariate_normal(mean=[0] * dims, cov=0.2 * np.eye(dims), size=Npts)
    return points + noise

max_stars_repo_name: rmsandu/Ellipsoid-F

max_issues_repo_name: rmsandu/Ellipsoid-F

max_forks_repo_name: raluca-san/Elipsoid-Fit

reason: import numpy,from scipy,import cvxpy


text:
import tensorflow as tf
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Conv2D, BatchNormalization, LeakyReLU
BatchNormalization._USE_V2_BEHAVIOR = False
from keras_contrib.layers import InstanceNormalization
from layers import ZeroPadding2D, ReflectionPadding2D, StridedConv

class Discriminator(Model):
    def __init__(self,
                 base_filters=32,
                 lrelu_alpha=0.2,
                 pad_type="reflect",
                 norm_type="batch"):
        super(Discriminator, self).__init__(name="Discriminator")
        if pad_type == "reflect":
            self.flat_pad = ReflectionPadding2D()
        elif pad_type == "constant":
            self.flat_pad = ZeroPadding2D()
        else:
            raise ValueError(f"pad_type not recognized {pad_type}")

        self.flat_conv = Conv2D(base_filters, 3)
        self.flat_lru = LeakyReLU(lrelu_alpha)
        self.strided_conv1 = StridedConv(base_filters * 2,
                                         lrelu_alpha,
                                         pad_type,
                                         norm_type)
        self.strided_conv2 = StridedConv(base_filters * 4,
                                         lrelu_alpha,
                                         pad_type,
                                         norm_type)
        self.conv2 = Conv2D(base_filters * 8, 3)

        if norm_type == "instance":
            self.norm = InstanceNormalization()
        elif norm_type == "batch":
            self.norm = BatchNormalization()

        self.lrelu = LeakyReLU(lrelu_alpha)

        self.final_conv = Conv2D(1, 3)

    def build(self, input_shape):
        super(Discriminator, self).build(input_shape)

    def call(self, x, training=False):
        x = self.flat_pad(x)
        x = self.flat_conv(x)
        x = self.flat_lru(x)
        x = self.strided_conv1(x, training=training)
        x = self.strided_conv2(x, training=training)
        x = self.conv2(x)
        x = self.norm(x, training=training)
        x = self.lrelu(x)
        x = self.final_conv(x)
        return x

if __name__ == "__main__":
    import numpy as np

    shape = (1, 256, 256, 3)
    nx = np.random.rand(*shape).astype(np.float32)
    t = tf.keras.Input(shape=nx.shape[1:], batch_size=nx.shape[0])
    tf.keras.backend.clear_session()
    sc = StridedConv(t.shape[-1])
    out = sc(t)
    sc.summary()
    print(f"Input  Shape: {t.shape}")
    print(f"Output Shape: {out.shape}")
    print("\n" * 2)

    tf.keras.backend.clear_session()
    d = Discriminator()
    out = d(t)
    d.summary()
    print(f"Input  Shape: {t.shape}")
    print(f"Output Shape: {out.shape}")


text:
import numpy as np
import graph_tool.all as gt
import scipy.sparse as sp
from sklearn.cluster import KMeans
from sklearn.preprocessing import MinMaxScaler

class NB_model:

    def __init__(self, g):
        self.g = g
        self.vals, self.vecs = self.calc_vertex_embedding()
        self.kmeans_labels, self.k = self.calc_kmeans()
        self.global_cen = None
        self.local_cen = None
        self.local_cen_com = None
        self.centrality = self.build_centrality()

    def calc_vertex_embedding(self):
        g = self.g

        bool_visited = g.new_vp("bool")
        bool_visited.a = g.vp.visited.a.astype(bool)
        g.set_vertex_filter(bool_visited, inverted=True)

        A = gt.adjacency(g)
        D = sp.diags([v.out_degree() for v in g.vertices()])
        Z = sp.csr_matrix((g.num_vertices(), g.num_vertices()))
        I = sp.eye(g.num_vertices(), g.num_vertices())
        up = sp.hstack((A, I - D))
        down = sp.hstack((I, Z))
        B = sp.vstack((up, down))

        if g.num_vertices() <= 1:
            return 1, np.zeros((B.shape))

        counter = 0
        while True:
            try:
                self.vals, self.vecs = sp.linalg.eigs(
                    B, k=min(B.shape[1] / 2, 6), which="LM", tol=0.01, maxiter=5000)
                break
            except:
                counter += 1
                if counter == 5:
                    break

        g.clear_filters()
        return self.vals, self.vecs

    def calc_kmeans(self):
        g = self.g
        vals = self.vals
        vecs = self.vecs
        c = 0.0
        k = 0
        if not np.iscomplex(vals[-1]):
            c = np.real(vals[-1] ** 0.5)
            for val in vals[::-1]:
                if not np.iscomplex(val) and np.real(val) > c:
                    k += 1
                else:
                    break
        else:
            k = 1

        k = min(k, 10)
        emb = np.real(vecs[:g.num_vertices(), -1:-1 - k - 1:-1])
        # community
        kmeans_labels = KMeans(init='k-means++', n_clusters=k,
                                    random_state=42).fit(emb).labels_
        g.vp.kmeans = g.new_vp("int", 0)
        g.vp.kmeans.a = kmeans_labels

        return kmeans_labels, k

    def calc_global_centrality(self):
        g = self.g
        vecs = self.vecs
        bool_visited = g.new_vp("bool")
        bool_visited.a = g.vp.visited.a.astype(bool)
        g.set_vertex_filter(bool_visited, inverted=True)

        global_cen = np.real(vecs[:g.num_vertices(), -1])
        self.global_cen = np.zeros((g.num_vertices(ignore_filter=True), 1))
        vertexMap = [int(v) for v in g.vertices()]
        for ind, val in enumerate(global_cen):
            self.global_cen[vertexMap[ind]] = abs(val)

        g.clear_filters()
        return self.global_cen

    # centrality for each community
    def calc_local_centrality(self):
        g = self.g
        k = self.k

        local_cen = g.new_vp("double", 0.0)
        local_cen_com = [g.new_vp("double", 0.0)
                         for i in range(0, k)]
        local_cen_a = local_cen.a
        local_cen_com_a = [cen.a for cen in local_cen_com]

        for i in range(0, k):
            free_nodes_com_p = g.new_vp("bool")
            free_nodes_com_p.a = np.logical_and(
                np.logical_not(g.vp.visited.a), g.vp.kmeans.a == i)
            g.set_vertex_filter(free_nodes_com_p)
            if g.num_vertices() <= 1:
                continue

            vertexMap = list(g.vertices())
            A = gt.adjacency(g)
            D = sp.diags([v.out_degree() for v in vertexMap])
            Z = sp.csr_matrix((g.num_vertices(), g.num_vertices()))
            I = sp.eye(g.num_vertices(), g.num_vertices())
            up = sp.hstack((A, I - D))
            down = sp.hstack((I, Z))
            B = sp.vstack((up, down))

            counter = 0
            while True:
                try:
                    vals, vecs = sp.linalg.eigs(B, k=1, which="LM", tol=0.005, maxiter=5000)
                    break
                except:
                    counter += 1
                    if counter == 5:
                        break

            for index, cenVal in enumerate(np.real(vecs[:g.num_vertices()])):
                local_cen_a[int(vertexMap[index])] = abs(cenVal)
                local_cen_com_a[i][int(vertexMap[index])] = abs(cenVal)
            g.clear_filters()

        self.local_cen = local_cen_a
        self.local_cen_com = local_cen_com
        return self.local_cen

    def build_centrality(self, is_rebuild=False):
        if is_rebuild:
            self.calc_vertex_embedding()

        global_centrality = self.calc_global_centrality()
        local_centrality = self.calc_local_centrality()
        local_centrality = np.array([[l] for l in local_centrality])
        state = np.hstack((local_centrality, global_centrality, self.kmeans_labels.reshape((self.g.num_vertices(ignore_filter=True), 1))))
        self.g.clear_filters()
        return MinMaxScaler(feature_range=(0, 1)).fit_transform(state)

# g = gt.collection.data["polbooks"]
# g.vp.visited = g.new_vp("bool", False)
# nbm = NB_model(g)
# pos = gt.sfdp_layout(g)
# cen = g.new_vp("double", 0.0)
# cen.a = nbm.global_cen.reshape((105))*30
# print cen.a
# gt.graph_draw(g, pos=pos, vertex_size=cen)

# for v in range(20):
#     g.vp.visited[g.vertex(v)] = True

# g.set_vertex_filter(g.vp.visited, inverted=True)

# nbm.calc_vertex_embedding()
# nbm.build_centrality()
# cen = g.new_vp("double", 0.0)
# cen.a = nbm.global_cen.reshape((105))*30
# g.set_vertex_filter(g.vp.visited, inverted=True)
# print cen.a
# gt.graph_draw(g, pos=pos, vertex_size=cen)

max_stars_repo_name: khurshedmemon/DQN-UN-TL

max_issues_repo_name: khurshedmemon/DQN-UN-TL

max_forks_repo_name: khurshedmemon/DQN-UN-TL

reason: import numpy,import scipy


text:
'''
Units tests for the cpptraj wrappers for nasqm
'''
import os
import pytest
import numpy as np
from pynasqm.nmr.nmrgroupgroup import NMRGroupGroup
import pynasqm.userinput as nasqm_user_input
import pynasqm.inputceon as inputceon

def setup_module(module):
    '''
    Switch to test directory
    '''
    os.chdir("tests/nmrTests/nmrGroupGroup")

def teardown_module(module):
    '''
    Return to main directory
    '''
    os.chdir("../../..")

def test_nmrgroupgroup_writer1():
    restricted_atoms1 = [[4055, 4056]]
    restricted_atoms2 = [[4065, 4066]]
    desired_distance = [4]
    writer = NMRGroupGroup(restricted_atoms1, restricted_atoms2, desired_distance)
    result_file = "nmr_groupgroup1.dist"
    answer_file = "nmr_groupgroup_test1.dist"
    writer.write_to(result_file)
    result = open(result_file, 'r').read()
    answer = open(answer_file, 'r').read()
    assert result == answer

def test_nmrgrougroup_writer2():
    restricted_atoms1 = [[4055, 4056], [4055, 4056]]
    restricted_atoms2 = [[4065, 4066], [5678, 5679]]
    desired_distance = [4,5]
    writer = NMRGroupGroup(restricted_atoms1, restricted_atoms2, desired_distance)
    result_file = "nmr_groupgroup2.dist"
    answer_file = "nmr_groupgroup_test2.dist"
    writer.write_to(result_file)
    result = open(result_file, 'r').read()
    answer = open(answer_file, 'r').read()
    assert result == answer

max_stars_repo_name: PotentialParadox/pynasqm

max_issues_repo_name: PotentialParadox/pynasqm

max_forks_repo_name: PotentialParadox/pynasqm


text:
import mgcpy.benchmarks.simulations as sims
import numpy as np
import pytest
from mgcpy.independence_tests.kendall_spearman import KendallSpearman

def test_kendall_spearman():
    # Against a randomly defined data set
    X = np.array([1.1728, 2.4941, 2.4101, 0.1814, 1.1978, 1.5806, 1.2504,
                  1.9706, 1.8839, 0.8760])[:, np.newaxis]
    Y = np.array([3.2311, 12.1113, 11.1350, 1.1989, 3.3127, 4.8580, 3.4917,
                  7.1748, 6.5792, 2.4012])[:, np.newaxis]
    kspear = KendallSpearman(None)
    kspear2 = KendallSpearman(None, 'spearman')

    test_stat1 = kspear.test_statistic(X, Y)[0]
    test_stat2 = kspear2.test_statistic(X, Y)[0]

    assert np.round(test_stat1, decimals=2) == 1.00
    assert np.round(test_stat2, decimals=2) == 1.00

    # Against linear simulations
    np.random.seed(0)
    X, Y = sims.linear_sim(100, 1)
    kspear = KendallSpearman(None)
    kspear2 = KendallSpearman(None, 'spearman')

    assert kspear.get_name() == 'kendall'
    assert kspear2.get_name() == 'spearman'

    test_stat1 = kspear.test_statistic(X, Y)[0]
    test_stat2 = kspear2.test_statistic(X, Y)[0]

    assert np.round(test_stat1, decimals=2) == 0.33
    assert np.round(test_stat2, decimals=2) == 0.48

max_stars_repo_name: neurodata-papers/mgc

max_issues_repo_name: neurodata-papers/mgc

max_forks_repo_name: neurodata-papers/mgc


text:
import getpass
from getpass import getuser
from typing import Callable, Tuple
import pdb
import random
from copy import deepcopy
import re

import numpy as np

# Get creation reactions for x mRNAs and y miRNAs, using the base model.
def base_creation_and_decay_reactions(ks, mus):
    string_list = ["#Reactions\n"]
    reaction_count = 0

    for RNA in range(len(ks)):
        string_list.append("R{}:\n".format(reaction_count))
        string_list.append("    $pool > r{}\n".format(RNA))
        string_list.append("    k{}\n\n".format(RNA))
        reaction_count += 1

    for RNA in range(len(mus)):
        string_list.append("R{}:\n".format(reaction_count))
        string_list.append("    r{} > $pool\n".format(RNA))
        string_list.append("    mu{} * r{}\n\n".format(RNA, RNA))
        reaction_count += 1

    string = ''.join(string_list)
    return string

# Get a string setting the creation and decay reactions for every
# RNA and every burst factor in a burst network.
# ks is the array of creation parameters.
# mus is the array of decay parameters,
def burst_creation_and_decay_reactions(n):
    string_list = ["#Reactions\n"]

    for RNA in range(n):
        # The reaction from the off burst state to the on burst state.
        # Reaction number
        string_list.append("R{}:\n".format(3*RNA))
        # Reaction equation
        string_list.append("    OFF{0} > ON{0}\n".format(RNA))
        # Reaction rate
        string_list.append("    alpha{0} * OFF{0}\n".format(RNA))

        # The reaction from the on burst state to the off burst state.
        # Reaction number
        string_list.append("R{}:\n".format(3*RNA+1))
        string_list.append("    ON{0} > OFF{0}\n".format(RNA))
        string_list.append("    beta{0} * ON{0}\n".format(RNA))

        # Creation of new RNA
        string_list.append("R{}:\n".format(3*RNA+2))
        string_list.append("    ON{0} > ON{0} + r{0}\n".format(RNA))
        string_list.append("    k{0} * ON{0}\n".format(RNA))

    # Degradation of RNA
    reaction_count = 3 * n
    for RNA in range(n):
        string_list.append("R{}:\n".format(RNA + reaction_count))
        string_list.append("    r{} > $pool\n".format(RNA))
        string_list.append("    mu{0} * r{0}\n\n".format(RNA))

    string = ''.join(string_list)
    return string

# Get a string setting all of the gamma interactions
# in the network.
# gammas is a dictionary containing the interactions and their strengths
# x is the number of mRNAs in the network
# reaction_count is the number of reactions that have already been added to the file
def get_gamma_reactions_string(x, gammas, reaction_count):
    string_list = []
    # For every mRNA
    for mRNA in gammas:
        # Not all values are guaranteed to be instantiated
        for miRNA, rate in gammas[mRNA].items():
            if rate > 0:
                string_list.append("R{}:\n".format(reaction_count))
                string_list.append("    r{} + r{} > $pool\n".format(mRNA, miRNA))
                string_list.append("    r{0} * r{1} * gammar{0}r{1}\n\n".format(mRNA, miRNA))
                reaction_count += 1
    string = ''.join(string_list)
    return string

# Get a string setting the initial species amounts
# for n RNAs, where n is the number of RNAs in the network.
def get_base_species_amounts(ks):
    string_list = ["#Species\n"]

    for RNA in range(len(ks)):
        # If the species has not been knocked out, set its initial value to 0.
        amount = 10 if ks[RNA] > 0 else 0
        string_list.append("    r{0} = {1}\n".format(RNA, amount))
    string = ''.join(string_list)
    return string

# Get a string setting the initial species amounts for
# n burst variables, where n is the number of RNAs in the network.
def get_burst_species_amounts(ks):
    string_list = []
    for RNA in range(len(ks)):
        # All species start in off state.
        string_list.append("    OFF{} = 1\n".format(RNA))
        string_list.append("    ON{} = 0\n".format(RNA))
    string_list.append("\n")
    string = ''.join(string_list)
    return string

# Get a string setting the creation and decay parameters for every RNA
# in the network.
# ks is an array giving the creation parameters.
# mus is an array giving the decay parameters.
def base_arrival_and_decay_parameters(ks: np.ndarray, mus: np.ndarray) -> str:
    string_list = ["#Parameters\n"]
    # parameter_string = "#Parameters\n"
    # For every RNA in th e network
    for RNA in range(len(ks)):
        string_list.append("    k{} = {}\n".format(RNA, ks[RNA]))
        string_list.append("    mu{} = {}\n".format(RNA, mus[RNA]))
    string = ''.join(string_list)
    return string

# Get a string setting all gamma interaction parameters for the network.
# x is the number of mRNAs in the network.
# gammas is a dictionary containing the gamma interactions.
def gamma_parameters(x, gammas):
    string_list = []
    # parameter_string = ""
    for RNA in gammas:
        for miRNA, rate in enumerate(gammas[RNA]):
            if rate > 0:
                string_list.append("    gammar{0}r{1} = {2}\n".format(RNA, miRNA, rate))
    string = ''.join(string_list)
    return string

# Get parameter strings for a set of burst modes
def burst_creation_and_decay_parameters(alphas, betas):
    string_list = []
    for species in range(len(alphas)):
        alpha = alphas[species]
        beta = betas[species]
        string_list.append("    alpha{0} = {1}".format(species, alpha))
        string_list.append("    beta{0} = {1}\n".format(species, beta))

    string = ''.join(string_list)
    return string

def get_base_network_writer(x: int, y: int):
    def file_writer(ks: np.ndarray, mus: np.ndarray, gammas: np.ndarray,
                    file_name: str):
        out_string_list = ["#Random base network with {0} mRNAs and {1} miRNAs\n"]

        # Remove any knockouts
        new_x, new_ks, new_mus, new_gammas = remove_knockouts(x, ks, mus, gammas)
        n = len(new_ks)

        # Reaction Strings
        out_string_list.append(base_creation_and_decay_reactions(new_ks, new_mus))
        out_string_list.append(get_gamma_reactions_string(new_x, new_gammas, 2 * n))

        # Species Strings
        out_string_list.append(get_base_species_amounts(ks))

        # Parameter Strings
        out_string_list.append(base_arrival_and_decay_parameters(new_ks, new_mus))
        out_string_list.append(gamma_parameters(new_x, new_gammas))

        out_string = ''.join(out_string_list)

        network_file = open("/home/" + getpass.getuser() + "/Stochpy/pscmodels/" + file_name + ".psc", "w")
        network_file.write(out_string)
        network_file.close()

    return file_writer

def get_burst_network_writer(x: int, y: int,
                             alphas: np.ndarray, betas: np.ndarray) -> Callable[..., type(None)]:

    def file_writer(ks: np.ndarray, mus: np.ndarray, gammas: np.ndarray, file_name: str) -> type(None):
        n = x + y
        out_string_list = ["#Random burst network with {0} mRNAs and {1} miRNAs".format(x, y)]
        # out_string = "#Random burst network of size " + str(x) + ", " + str(y) + ".\n"

        # Reaction Strings
        out_string_list.append(burst_creation_and_decay_reactions(n))
        # out_string += burst_creation_and_decay_reactions(n)
        # out_string += get_gamma_reactions_string(x, gammas, 4 * (x + y))
        out_string_list.append(get_gamma_reactions_string(x, gammas, 4 * (x + y)))

        # Species Strings
        # out_string += get_base_species_amounts(ks)
        out_string_list.append(get_base_species_amounts(ks))
        # out_string += get_burst_species_amounts(ks)
        out_string_list.append(get_burst_species_amounts(ks))

        # Parameter Strings
        # out_string += base_arrival_and_decay_parameters(ks, mus)
        out_string_list.append(base_arrival_and_decay_parameters(ks, mus))
        # out_string += gamma_parameters(x, gammas)
        out_string_list.append(gamma_parameters(x, gammas))
        # out_string += burst_creation_and_decay_parameters(alphas, betas)
        out_string_list.append(burst_creation_and_decay_parameters(alphas, betas))

        string = ''.join(out_string_list)
        # assert string == out_string
        network_file = open("/home/" + getpass.getuser() + "/Stochpy/pscmodels/" + file_name + ".psc", "w")
        network_file.write(string)
        network_file.close()

    return file_writer

def get_complex_network_writer(x: int, y: int,
                               alphas: np.ndarray, betas: np.ndarray, deltas: np.ndarray) -> Callable[..., type(None)]:

    def file_writer(ks: np.ndarray, mus: np.ndarray, gammas: np.ndarray, file_name: str) -> type(None):
        n = x + y

        out_string_list = ["#Random burst network with {0} mRNAs and {1} miRNAs".format(x, y)]
        # out_string = "#Random burst network of size " + str(x) + ", " + str(y) + ".\n"

        # Reaction Strings
        out_string_list.append(burst_creation_and_decay_reactions(n))
        # out_string += burst_creation_and_decay_reactions(n)
        # out_string += get_gamma_reactions_string(x, gammas, 4 * (x + y))
        out_string_list.append(get_gamma_reactions_string(x, gammas, 4 * (x + y)))

        # Species Strings
        # out_string += get_base_species_amounts(ks)
        out_string_list.append(get_base_species_amounts(ks))
        # out_string += get_burst_species_amounts(ks)
        out_string_list.append(get_burst_species_amounts(ks))

        # Parameter Strings
        # out_string += base_arrival_and_decay_parameters(ks, mus)
        out_string_list.append(base_arrival_and_decay_parameters(ks, mus))
        # out_string += gamma_parameters(x, gammas)
        out_string_list.append(gamma_parameters(x, gammas))
        # out_string += burst_creation_and_decay_parameters(alphas, betas)
        out_string_list.append(burst_creation_and_decay_parameters(alphas, betas))

        string = ''.join(out_string_list)
        # assert string == out_string
        network_file = open("/home/" + getpass.getuser() + "/Stochpy/pscmodels/" + file_name + ".psc", "w")
        network_file.write(string)
        network_file.close()

    return file_writer

# Knockout species (k_i = 0) must be genuinely removed, or the simulations
# won't run as well.
def remove_knockouts(x: int, ks: np.ndarray, mus: np.ndarray, gammas: dict,
                     alphas: np.ndarray = None, betas: np.ndarray = None, deltas: dict = None) -> tuple:
    new_x = 0
    new_ks = []
    new_mus = []
    new_gammas = {}
    if alphas:
        new_alphas = []
        new_betas = []
        if deltas:
            new_deltas = []
    non_knockouts = np.nonzero(ks)[0]
    for species in non_knockouts:
        # Keep track of how many mRNAs are in the new network.
        if species < x:
            new_x += 1
        new_ks.append(ks[species])
        new_mus.append(mus[species])
        new_gammas[species] = {x: gammas[species][x] for x in non_knockouts}

    new_ks = np.array(new_ks)
    new_mus = np.array(new_mus)
    return new_x, new_ks, new_mus, new_gammas

max_stars_repo_name: VincentLuczkow/cerna

max_issues_repo_name: VincentLuczkow/cerna

max_forks_repo_name: vluzko/cerna