text:
"""Module implementing point transformations and their matrices."""
import numpy as np

def axis_angle_rotation(axis, angle, point=None, deg=True):
    r"""Return a 4x4 matrix for rotation about any axis by given angle.

    Rotations around an axis that contains the origin can easily be
    computed using Rodrigues' rotation formula. The key quantity is
    the ``K`` cross product matrix for the unit vector ``n`` defining
    the axis of the rotation:

             /   0  -nz   ny \
        K =  |  nz    0  -nx |
             \ -ny   nx    0 /

    For a rotation angle ``phi`` around the vector ``n`` the rotation
    matrix is given by

        R = I + sin(phi) K  + (1 - cos(phi)) K^2

    where ``I`` is the 3-by-3 unit matrix and ``K^2`` denotes the matrix
    square of ``K``.

    If the rotation axis doesn't contain the origin, we have to first
    shift real space to transform the axis' ``p0`` reference point into
    the origin, then shift the points back after rotation:

        p' = R @ (p - p0) + p0 = R @ p + (p0 - R @ p0)

    This means that the rotation in general consists of a 3-by-3
    rotation matrix ``R``, and a translation given by
    ``b = p0 - R @ p0``. These can be encoded in a 4-by-4 transformation
    matrix by filling the 3-by-3 leading principal submatrix with ``R``,
    and filling the top 3 values in the last column with ``b``.

    Parameters
    ----------
    axis : 3-length sequence
        The direction vector of the rotation axis. It need not be a
        unit vector, but it must not be a zero vector.

    angle : float
        Angle of rotation around the axis. The angle is defined as a
        counterclockwise rotation when facing the normal vector of the
        rotation axis. Passed either in degrees or radians depending on
        the value of ``deg``.

    point : 3-length sequence, optional
        The origin of the rotation (a reference point through which the
        rotation axis passes). By default the rotation axis contains the
        origin.

    deg : bool, optional
        Whether the angle is specified in degrees. ``False`` implies
        radians.

    Examples
    --------
    Generate a transformation matrix for rotation around a cube's body
    diagonal by 120 degrees.

    >>> import numpy as np
    >>> from pyvista import transformations
    >>> trans = transformations.axis_angle_rotation([1, 1, 1], 120)

    Check that the transformation cycles the cube's three corners.

    >>> corners = np.array([
    ...     [1, 0, 0],
    ...     [0, 1, 0],
    ...     [0, 0, 1],
    ... ])
    >>> rotated = transformations.apply_transformation_to_points(trans, corners)
    >>> np.allclose(rotated, corners[[1, 2, 0], :])
    True

    """
    if deg:
        # convert to radians
        angle *= np.pi / 180

    # return early for no rotation; play it safe and check only exact equality
    if angle % (2 * np.pi) == 0:
        return np.eye(4)

    axis = np.asarray(axis, dtype='float64')
    if axis.shape != (3,):
        raise ValueError('Axis must be a 3-length array-like.')
    if point is not None:
        point = np.asarray(point)
        if point.shape != (3,):
            raise ValueError('Rotation center must be a 3-length array-like.')

    # check and normalize
    axis_norm = np.linalg.norm(axis)
    if np.isclose(axis_norm, 0):
        raise ValueError('Cannot rotate around zero vector axis.')
    if not np.isclose(axis_norm, 1):
        axis = axis / axis_norm

    # build Rodrigues' rotation matrix
    K = np.zeros((3, 3))
    K[[2, 0, 1], [1, 2, 0]] = axis
    K += -K.T
    R = np.eye(3) + np.sin(angle) * K + (1 - np.cos(angle)) * K @ K
    augmented = np.eye(4)
    augmented[:-1, :-1] = R

    if point is not None:
        # rotation of point p would be R @ (p - point) + point
        # which is R @ p + (point - R @ point)
        augmented[:-1, -1] = point - R @ point

    return augmented

def reflection(normal, point=None):
    """Return a 4x4 matrix for reflection across a normal about a point.

    Projection to a unit vector ``n`` can be computed using the dyadic
    product (or outer product) ``P`` of ``n`` with itself, which is a
    3-by-3 symmetric matrix.

    Reflection across a plane that contains the origin amounts to
    reversing the components of real space points that are perpendicular
    to the reflection plane. This gives us the transformation ``R``
    acting on a point ``p`` as

        p' = R @ p = p - 2 P @ p = (I - 2 P) @ p

    so the reflection's transformation matrix is the unit matrix minus
    twice the dyadic product ``P``.

    If additionally we want to compute a reflection to a plane that does
    not contain the origin, we can we can first shift every point in
    real space by ``-p0`` (if ``p0`` is a point that lies on the plane)

        p' = R @ (p - p0) + p0 = R @ p + (p0 - R @ p0)

    This means that the reflection in general consists of a 3-by-3
    reflection matrix ``R``, and a translation given by
    ``b = p0 - R @ p0``. These can be encoded in a 4-by-4 transformation
    matrix by filling the 3-by-3 leading principal submatrix with ``R``,
    and filling the top 3 values in the last column with ``b``.

    Parameters
    ----------
    normal : 3-length sequence
        The normal vector of the reflection plane. It need not be a unit
        vector, but it must not be a zero vector.

    point : 3-length sequence, optional
        The origin of the reflection (a reference point through which
        the reflection plane passes). By default the reflection plane
        contains the origin.

    Examples
    --------
    Generate a transformation matrix for reflection over the XZ plane.

    >>> import numpy as np
    >>> from pyvista import transformations
    >>> trans = transformations.reflection([0, 1, 0])

    Check that the reflection transforms corners of a cube among one
    another.

    >>> verts = np.array([
    ...     [ 1, -1,  1],
    ...     [-1, -1,  1],
    ...     [-1, -1, -1],
    ...     [-1, -1,  1],
    ...     [ 1,  1,  1],
    ...     [-1,  1,  1],
    ...     [-1,  1, -1],
    ...     [-1,  1,  1],
    ... ])
    >>> mirrored = transformations.apply_transformation_to_points(trans, verts)
    >>> np.allclose(mirrored, verts[[np.r_[4:8, 0:4]], :])
    True

    """
    normal = np.asarray(normal, dtype='float64')
    if normal.shape != (3,):
        raise ValueError('Normal must be a 3-length array-like.')
    if point is not None:
        point = np.asarray(point)
        if point.shape != (3,):
            raise ValueError('Plane reference point must be a 3-length array-like.')

    # check and normalize
    normal_norm = np.linalg.norm(normal)
    if np.isclose(normal_norm, 0):
        raise ValueError('Plane normal cannot be zero.')
    if not np.isclose(normal_norm, 1):
        normal = normal / normal_norm

    # build reflection matrix
    projection = np.outer(normal, normal)
    R = np.eye(3) - 2 * projection
    augmented = np.eye(4)
    augmented[:-1, :-1] = R

    if point is not None:
        # reflection of point p would be R @ (p - point) + point
        # which is R @ p + (point - R @ point)
        augmented[:-1, -1] = point - R @ point

    return augmented

def apply_transformation_to_points(transformation, points, inplace=False):
    """Apply a given transformation matrix (3x3 or 4x4) to a set of points.

    Parameters
    ----------
    transformation : np.ndarray
        Transformation matrix of shape (3, 3) or (4, 4).

    points : np.ndarray
        Array of points to be transformed of shape (N, 3).

    inplace : bool, optional
        Updates points in-place while returning nothing.

    Returns
    -------
    numpy.ndarray
        Transformed points.

    Examples
    --------
    Scale a set of points in-place.

    >>> import numpy as np
    >>> import pyvista
    >>> from pyvista import examples
    >>> points = examples.load_airplane().points
    >>> points_orig = points.copy()
    >>> scale_factor = 2
    >>> tf = scale_factor * np.eye(4)
    >>> tf[3, 3,] = 1
    >>> pyvista.transformations.apply_transformation_to_points(tf, points, inplace=True)
    >>> assert np.all(np.isclose(points, scale_factor * points_orig))

    """
    transformation_shape = transformation.shape
    if transformation_shape not in ((3, 3), (4, 4)):
        raise ValueError('`transformation` must be of shape (3, 3) or (4, 4).')

    if points.shape[1] != 3:
        raise ValueError('`points` must be of shape (N, 3).')

    if transformation_shape[0] == 4:
        # Divide by scale factor when homogeneous
        transformation /= transformation[3, 3]

        # Add the homogeneous coordinate
        # `points_2` is a copy of the data, not a view
        points_2 = np.empty((len(points), 4))
        points_2[:, :-1] = points
        points_2[:, -1] = 1
    else:
        points_2 = points

    # Paged matrix multiplication. For arrays with ndim > 2, matmul assumes
    # that the matrices to be multiplied lie in the last two dimensions.
    points_2 = (transformation[np.newaxis, :, :] @ points_2.T)[0, :3, :].T

    # If inplace, set the points
    if inplace:
        points[:] = points_2
    else:
        # otherwise return the new points
        return points_2

max_stars_repo_name: eino/pyvista

max_issues_repo_name: eino/pyvista

max_forks_repo_name: eino/pyvista


text:
#  Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
#  See https://llvm.org/LICENSE.txt for license information.
#  SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

import re
from typing import Iterable
import numpy as np

from _npcomp.mlir import ir

from npcomp.dialect import Numpy
from npcomp.exporter import *
from npcomp.types import *
from npcomp.tracing.context import *
from npcomp.tracing.emitters import *

class ModuleBuilder:
  """Builds an MLIR module by tracing functions."""

  def __init__(self, mlir_context=None, emitter_registry=None):
    self.context = mlir_context if mlir_context else ir.MLIRContext()
    self.module = self.context.new_module()
    self.helper = Numpy.DialectHelper(self.context, ir.OpBuilder(self.context))
    self.emitters = (emitter_registry
                     if emitter_registry else EmitterRegistry.create_default())

  def trace(self, *export_py_funcs: ExportPyFunction):
    """Traces exported py functions."""
    for export_py_func in export_py_funcs:
      assert isinstance(export_py_func, ExportPyFunction), (
          "Expected an exported python function (from the Exporter class)")
      tracer = FunctionTracer(self, export_py_func)
      with tracer:
        tracer.trace()

class FunctionTracer(TraceContext):
  """A trace of a single function."""
  __slots__ = [
      "module_builder",
      "epf",
      "_args_array_params",
      "_f",
      "_f_types",
      "_helper",
      "_mlir_m",
      "_mlir_c",
      "_python_args",
      "_result_array_params",
      "_traced_arrays",
      "_external_arrays",
  ]

  def __init__(self, module_builder: ModuleBuilder, epf: ExportPyFunction):
    super().__init__(desc="[trace of %s]" % epf.__name__)
    self.module_builder = module_builder
    self.epf = epf
    self._traced_arrays = {}  # Mapping of TracedArray to current consumer value
    self._external_arrays = {}  # Mapping of id to (ndarray, ir.Value)
    self._validate()

    # Alias some parent members for convenience.
    self._mlir_m = module_builder.module
    self._mlir_c = module_builder.context
    self._helper = module_builder.helper

    # Extract ArrayParams for all args and results.
    self._args_array_params = [
        ArrayParams.from_constraints(arg.constraints)
        for arg in self.epf.sig.args
    ]
    self._python_args = [None] * len(self._args_array_params)
    self._result_array_params = ArrayParams.from_constraints(
        self.epf.sig.result.constraints)

    # Create the MLIR function.
    self._f, self._f_types = self._create_mlir_function()
    self._create_trace_roots()

  def trace(self):
    # Invoke the python function with placeholders.
    # TODO: More sophisticated signature merging
    # TODO: Multiple results
    # TODO: Error reporting
    h = self._helper
    py_results = (self.epf.pyfunc(*self._python_args),)
    if len(py_results) != len(self._f_types):
      raise TracingError("Traced function returned != %d results: %r" % (
          len(self._f_types),
          py_results,
      ))

    # Narrow all results to the declared return types.
    return_operands = []
    for py_result, mlir_result_type in zip(py_results, self._f_types):
      mlir_result = self.get_traced_array_value(py_result)
      # narrow to declared result type.
      return_operands.extend(
          h.numpy_narrow_op(mlir_result_type, mlir_result).results)
    h.return_op(return_operands)

  def set_traced_array(self, traced_array, value):
    """Sets the current SSA value for a traced_array."""
    assert isinstance(traced_array, TracedArray)
    self._traced_arrays[traced_array] = value

  def get_traced_array_value(self, traced_array):
    if not isinstance(traced_array, TracedArray):
      # Generic import of external value. For now, we just treat these as
      # local consts.
      return self._get_external_array_value(traced_array)

    traced_value = self._traced_arrays.get(traced_array)
    if traced_value is None:
      raise TracingError("Unregistered traced array: %r", (traced_array,))
    return traced_value

  def _get_external_array_value(self, external_array):
    h = self._helper
    if not isinstance(external_array, np.ndarray):
      raise TracingError("Expected ndarray but got: %r" % (external_array,))
    found_it = self._external_arrays.get(id(external_array))
    if found_it:
      return found_it[1]
    # Import it.
    dense_attr = h.context.dense_elements_attr(external_array)
    const_value = h.constant_op(dense_attr.type, dense_attr).result
    self._external_arrays[id(external_array)] = (external_array, const_value)
    return const_value

  def _validate(self):
    if not all(
        arg.type_class == TypeClass.NdArray for arg in self.epf.sig.args):
      raise NotImplementedError("Non NdArray args: %r" % (self.epf.sig.args,))
    if not self.epf.sig.result.type_class == TypeClass.NdArray:
      raise NotImplementedError("Non NdArray result: %r" %
                                (self.epf.sig.result,))

  def _create_mlir_function(self):
    mlir_c = self._mlir_c
    mlir_m = self._mlir_m
    h = self._helper
    epf = self.epf
    f_args = [
        mlir_c.parse_type(ap.mlir_tensor_type_asm)
        for ap in self._args_array_params
    ]
    f_types = [
        mlir_c.parse_type(self._result_array_params.mlir_tensor_type_asm)
    ]
    h.builder.insert_before_terminator(mlir_m.first_block)
    f_type = h.function_type(f_args, f_types)
    f = h.func_op(epf.__name__, f_type, create_entry_block=True)
    return f, f_types

  def _create_trace_roots(self):
    entry_block = self._f.first_block
    for index, ap in enumerate(self._args_array_params):
      if ap is not None:
        ta = TracedArray(self)
        self.set_traced_array(ta, entry_block.args[index])
        self._python_args[index] = ta

  def _resolve_input_ssa_values(self, trace_values: Iterable[TraceValue]):
    """Resolves input python values to SSA values."""
    ssa_values = []
    for tv in trace_values:
      assert tv.type == TraceValueType.NDARRAY, (
          "Unsupported TraceValueType: %r" % tv.type)
      ssa_value = self.get_traced_array_value(tv.value)
      ssa_values.append(ssa_value)
    return ssa_values

  def _resolve_result_py_values(self,
                                trace_value_types: Iterable[TraceValueType],
                                ssa_values):
    """Resolves result SSA values to runtime python values."""
    assert len(trace_value_types) == len(ssa_values), (
        "Mismatched emitter declared result types and results")
    py_values = []
    for trace_value_type, ssa_value in zip(trace_value_types, ssa_values):
      assert trace_value_type == TraceValueType.NDARRAY, (
          "Unsupported TraceValueType: %r" % trace_value_type)
      py_value = TracedArray(self)
      self.set_traced_array(py_value, ssa_value)
      py_values.append(py_value)
    return py_values

  def _emit_invocation(self, emitter: FuncEmitter, invocation: TraceInvocation):
    tv_map = emitter.map_invocation(invocation)
    input_ssa_values = self._resolve_input_ssa_values(tv_map.input_trace_values)
    request = EmissionRequest(input_ssa_values,
                              dialect_helper=self._helper,
                              extra=tv_map.extra)
    result_ssa_values = emitter.emit(request)
    py_values = self._resolve_result_py_values(tv_map.result_trace_value_types,
                                               result_ssa_values)
    return emitter.map_results(py_values, tv_map.extra)

  def _handle_ufunc(self, ufunc, method, inputs, kwargs):
    emitter = self.module_builder.emitters.lookup_ufunc(ufunc, method)
    if not emitter:
      return NotImplemented
    invocation = TraceInvocation(inputs, kwargs, Protocol.UFUNC, method)
    return self._emit_invocation(emitter, invocation)

  def _handle_array_func(self, func, types, inputs, kwargs):
    emitter = self.module_builder.emitters.lookup_array_func(func)
    if not emitter:
      return NotImplemented
    invocation = TraceInvocation(inputs, kwargs, Protocol.ARRAY_FUNC)
    return self._emit_invocation(emitter, invocation)

  def _emit_slice_value(self, slice_element):
    h = self._helper
    if slice_element == None:
      return h.basicpy_singleton_op(h.basicpy_NoneType).result
    elif slice_element == Ellipsis:
      return h.basicpy_singleton_op(h.basicpy_EllipsisType).result
    elif isinstance(slice_element, int):
      return h.constant_op(h.index_type,
                           h.context.index_attr(slice_element)).result
    elif isinstance(slice_element, slice):
      return self._emit_slice_object(slice_element)
    else:
      # Assume array convertible.
      raise NotImplementedError(
          "TODO: Slicing with generic arrays not yet implemented")

  def _emit_slice_object(self, slice_object: slice):
    h = self._helper

    def emit_index(index):
      if index is None:
        return h.basicpy_singleton_op(h.basicpy_NoneType).result
      else:
        return h.constant_op(h.index_type,
                             h.context.index_attr(int(index))).result

    start = emit_index(slice_object.start)
    stop = emit_index(slice_object.stop)
    step = emit_index(slice_object.step)
    return h.basicpy_slot_object_make_op("slice", start, stop, step).result

  def _handle_array_getitem(self, array, key):
    h = self._helper
    array_value = self.get_traced_array_value(array)
    # Array slicing is always based on a tuple.
    slice_tuple = key if isinstance(key, tuple) else (key,)
    # Resolve and emit each slice element.
    slice_values = [self._emit_slice_value(elt) for elt in slice_tuple]
    result_value = h.numpy_get_slice_op(h.unknown_array_type, array_value,
                                        *slice_values).result
    result_array = TracedArray(self)
    self.set_traced_array(result_array, result_value)
    return result_array

if __name__ == "__main__":
  import doctest
  doctest.testmod()

ext: py

max_stars_repo_name: marbre/mlir-npcomp

max_issues_repo_name: marbre/mlir-npcomp

max_forks_repo_name: marbre/mlir-npcomp

reason: import numpy


text:
"""
    core/__audio.py

    helpful with development of programs
    that work with audio

    author: @alexzander
"""

# python
import os
import threading
import subprocess

# 3rd party
import pyttsx3 # pip install pyttsx3
# import sounddevice # pip install sounddevice
from gtts import gTTS # pip install gtts
from pydub import AudioSegment # pip install pydub
from playsound import playsound # pip install playsound
import speech_recognition as sr # pip install SpeechRecognition
from scipy.io.wavfile import write # pip install scipy

# core package (pip install python-core)
from core.system import *
from core.path__ import *
from core.aesthetics import *

def playaudio(path, __threaded=False):
    """
        plays sound threaded or not threaded
    """
    if type(path) != str:
        raise TypeError("param @absolute_path should be str.")
    if not os.path.isfile(path):
        raise ValueError(f"param {path} is not an absolute path.")

    if __threaded:
        threading.Thread(target=playsound, args=(path, )).start()
    else:
        playsound(path)

def save_text_to_speech(
    text: str,
    destination: str,
    language="en",
    slowmotion=False,
    __print=False
):
    if type(text) != str:
        text = str(text)
    if text == "":
        raise ValueError("param @text is void-string.")

    # @dst
    if type(destination) != str:
        raise TypeError("param @dst should be type str.")
    if not destination.endswith(".mp3") and not destination.endswith(".wav"):
        raise ValueError("param @dst should end with .mp3")
    if not is_abs(destination):
        raise ValueError

    # @lang
    if type(language) != str:
        language = str(language)

    speech = gTTS(text, lang=language, slow=slowmotion)
    speech.save(destination)

    if __print:
        message = "Text: '{}'\nwas saved to speech here: {}\n{}".format(yellow(text), cyan(destination), green_bold("successfully."))
        print(message)

def save_recorded_mic(destination, duration):
    # validation
    if type(destination) != str:
        raise TypeError("param @dst should be type str.")
    if type(duration) != int:
        raise TypeError("param @duration should be type int.")

    if not is_abs(destination):
        raise ValueError("param @dst should be absolute path.")

    if not destination.endswith(".wav"):
        raise ValueError("param @dst should end with .wav extension.")

    if duration <= 0:
        raise ValueError("param @duration should be bigger than 0.")
    # /validation

    try:
        sample_rate = 44100
        print("you have {} seconds to speak:".format(duration))
        print("recording...")

        myrecording = sounddevice.rec(int(duration * sample_rate), samplerate=sample_rate, channels=2)
        sounddevice.wait()

        write(destination, sample_rate, myrecording)

        print("recording saved to path: {}.".format(destination))

    except OSError as error:
        print(type(error))
        print(error)
        print("turn OFF your anti-virus application.")

def speech_to_text_from_mic():
    """ listens to your voice, if understands:
        return speech to text else return None.
    """
    recognizer = sr.Recognizer()
    microphone = sr.Microphone()

    print("listening...")
    with microphone as audio:
        recognizer.adjust_for_ambient_noise(audio)
        recorded_audio = recognizer.listen(audio)
    try:
        text = recognizer.recognize_google(recorded_audio)
    except sr.UnknownValueError:
        print("your speech was incomprehensible. None was returned.")
        return None
    return text

def speech_to_text_from_file(path):
    if not is_file(path):
        raise ValueError
    if not path.endswith(".mp3") and not path.endswith(".wav"):
        raise ValueError

    r = sr.Recognizer()
    with sr.AudioFile(path) as audio_file:
        data = r.record(audio_file)
        try:
            text = r.recognize_google(data)
            return text
        except sr.UnknownValueError as error:
            print(error)
            print(error.message)
            print(type(error))

def speak_engine(message, speech_speed=150):
    if type(message) != str:
        raise TypeError("param @message should be type str.")
    if message == "":
        raise ValueError("param @message shouldn't be void-string.")

    engine = pyttsx3.init("sapi5")
    engine.setProperty("rate", speech_speed)
    engine.setProperty("volume", 0.8)
    voices = engine.getProperty("voices")
    engine.setProperty("voice", voices[1].id)
    engine.say(message)
    engine.runAndWait()

def modify_volume(src, dst, units=15):
    """ applicable only to .wav and .mp3 files. """
    # validation
    if type(src) != str:
        message = ConsoleColored("param @absolute_path should be type str.", "red", bold=1)
        raise TypeError(message)
        del message
    if type(units) != int:
        message = ConsoleColored("param @quantity should be type int.", "red", bold=1)
        raise TypeError(message)
        del message

    if not os.path.isfile(src):
        message = ConsoleColored("param @absolute_path is not an absolut path.", "red", bold=1)
        raise ValueError(message)
        del message

    if not src.endswith(".wav") and not src.endswith("mp3"):
        message = ConsoleColored("param @absolute_path should end with .wav or .mp3.", "red", bold=1)
        raise ValueError(message)
        del message
    # /validation

    extension = get_file_extension(src)
    if extension == "wav":
        audio_file = AudioSegment.from_wav(src)
    elif extension == "mp3":
        audio_file = AudioSegment.from_mp3(src)
    else:
        raise InvalidExtensionError("{} has {} extension".format(src, extension))

    filename = get_file_name(src)
    audio_file += units
    if units < 0:
        modified_file = dst + "\\{}_minus_{}.{}".format(filename, -units, extension)
    else:
        modified_file = dst + "\\{}_plus_{}.{}".format(filename, units, extension)

    audio_file.export(modified_file, extension)

    # code=0
    return True

# WE HAVE A PROBLEM WITH THESE
# WE CANT MAKE EXECUTABLE WITH THESE MODULES
# IT CAUSES SOME ERROR AND RECURSION ERROR
# LIBROSA HAS A PROBLEM TOO (its not finding a .txt in )

# watch for this D:\Alexzander__\programming\python\202020_order\error.png

# =======================================================
# import numpy as np
# import pandas as pd
# import matplotlib.pyplot as plt
# from glob import glob

# pip install librosa as administrator
# import librosa as lr

# def vizualize_audio(wav_file):
#     audio, sfreq = lr.load(wav_file)
#     __time = np.arange(0, len(audio)) / sfreq

#     fig, ax = plt.subplots()
#     ax.plot(__time, audio)
#     f_name = get_filename_plus_extension(wav_file)
#     ax.set(xlabel="Time in seconds", ylabel="Sound amplitude", title=f_name)
#     print("visualization of: {}\nwas {} in a gui window.".format(
#         blue_bold(f_name),
#         yellow_bold("opened")
#     ))
#     plt.show()
# =======================================================

import wave
import pyaudio
import time

def convert_file_to_wav(mp3_file: str, __print=False):
    _sep = get_path_sep(mp3_file)
    _items = mp3_file.split(_sep)
    destination_folder = _sep.join(_items[:-1])
    filename = _items[-1]
    full_dest = "{}/{}__converted_to_wav.wav".format(destination_folder, filename)

    try:
        # actual convertion
        # you need to have ffmpeg.exe file
        # in the environment variables to do this
        subprocess.call(["ffmpeg", "-i", mp3_file, full_dest])
        if __print:
            print(f"File: {mp3_file}")
            print_green("converted to WAV successfully")
            print(f"Location: {full_dest}")

    except BaseException as err:
        if __print:
            print("\nsome unexpected error:")
            print(err)
            print(type(TypeError))

def play_wav(wav_file):
    chunk = 1024
    wf = wave.open(wav_file, 'rb')
    p = pyaudio.PyAudio()

    stream = p.open(
        format = p.get_format_from_width(wf.getsampwidth()),
        channels = wf.getnchannels(),
        rate = wf.getframerate(),
        output = True
    )

    try:
        print("Currently playing:")
        print(wav_file)

        data = wf.readframes(chunk)
        while data:
            # writing to the stream is what *actually* plays the sound.
            stream.write(data)
            data = wf.readframes(chunk)

    except KeyboardInterrupt:
        print("music stopped from [ KeyboardInterrupt ]")
        pass

    stream.close()
    wf.close()
    p.terminate()

# TESTING
if __name__ == '__main__':
    # song = r"D:\Alexzander__\audio\best-songs\mrbot.mp3"
    # convert_file_to_wav(song, True)
    pass

reason: import numpy,from scipy


text:
"""
Created on Aug 15, 2011

@author: guillaume
"""

import scipy as sc
from scipy.linalg import expm2 as expm

from chemex.experiments.misc import correct_chemical_shift
from chemex.caching import lru_cache
from .liouvillian import compute_nz_eq, compute_base_liouvillians, compute_free_liouvillian, get_nz

@lru_cache()
def make_calc_observable(time_t1=0.0, b1_offset=0.0, b1_frq=0.0, b1_inh=0.0, b1_inh_res=5, carrier=0.0, ppm_to_rads=0.0,
                         _id=None):
    """
    Factory to make "calc_observable" function to calculate the intensity in presence
    of exchange after a CEST block.

    Parameters
    ----------
    time_t1 : float
        Duration of the CW block.
    b1_offset : float
        Frequency offset of the applied B1 field in Hz.
    b1_frq : float
        Strength of the applied B1 field in Hz.
    b1_inh : float
        B1 field inhomogeneity in Hz.
    b1_inh_res : int
        Resolution to model B1 field inhomogeneity.
    carrier : float
        Carrier position in rad/s.
    ppm_to_rads : float
        Conversion factor from ppm to rad/s
    id : tuple
        Some type of identification for caching optimization

    Returns
    -------
    out : function
        Calculate intensity after the CEST block

    """

    base_liouvillians, weights = compute_base_liouvillians(b1_offset, b1_frq, b1_inh, b1_inh_res)

    @lru_cache(5)
    def _calc_observable(pb=0.0, kex=0.0, dw=0.0, r_nz=1.5, r_nxy=0.0, dr_nxy=0.0, cs=0.0):
        """
        Calculate the intensity in presence of exchange after a CEST block assuming
        initial intensity of 1.0.

        Parameters
        ----------
        pb : float
            Fractional population of state B,
            0.0 for 0%, 1.0 for 100%
        kex : float
            Exchange rate between state A and B in /s.
        dw : float
            Chemical shift difference between states A and B in rad/s.
        r_nz : float
            Longitudinal relaxation rate of state {a,b} in /s.
        r_nxy : float
            Transverse relaxation rate of state a in /s.
        dr_nxy : float
            Transverse relaxation rate difference between states a and b in /s.
        cs : float
            Resonance position in rad/s.

        Returns
        -------
        out : float
            Intensity after the CEST block

        """

        if abs(b1_offset) >= 10000.0:

            return 1.0 - pb

        else:

            dw *= ppm_to_rads
            mag_eq = compute_nz_eq(pb)
            exchange_induced_shift, _ = correct_chemical_shift(pb=pb, kex=kex, dw=dw,
                                                               r_ixy=r_nxy, dr_ixy=dr_nxy)
            wg = (cs - carrier) * ppm_to_rads - exchange_induced_shift

            liouvillians = base_liouvillians + compute_free_liouvillian(pb=pb, kex=kex, dw=dw,
                                                                        r_nxy=r_nxy, dr_nxy=dr_nxy,
                                                                        r_nz=r_nz, cs_offset=wg)

            propagator = sc.zeros_like(liouvillians[0])
            for liouvillian, weight in zip(liouvillians, weights):
                propagator += weight * expm(liouvillian * time_t1)
            propagator /= sum(weights)

            magz_a, _ = get_nz(sc.dot(propagator, mag_eq))

            return magz_a

    def calc_observable(i0=0.0, **kwargs):
        """
        Calculate the intensity in presence of exchange after a CEST block.

        Parameters
        ----------
        i0 : float
            Initial intensity.

        Returns
        -------
        out : float
            Intensity after the CEST block

        """

        return i0 * _calc_observable(**kwargs)

    return calc_observable

reason: import scipy,from scipy


text:
"""
.. codeauthor:: David Zwicker <david.zwicker@ds.mpg.de>
"""

import numpy as np
import pytest

from pde import CahnHilliardPDE, CartesianGrid, DiffusionPDE, ScalarField, UnitGrid
from pde.tools.misc import skipUnlessModule

from droplets import LengthScaleTracker, SphericalDroplet
from droplets.emulsions import EmulsionTimeCourse

@skipUnlessModule("h5py")
def test_emulsion_tracker(tmp_path):
    """test using the emulsions tracker"""
    path = tmp_path / "test_emulsion_tracker.hdf5"

    d = SphericalDroplet([4, 4], 3)
    c = d.get_phase_field(UnitGrid([8, 8]))

    pde = CahnHilliardPDE()

    e1 = EmulsionTimeCourse()
    tracker = e1.tracker(filename=path)
    pde.solve(c, t_range=1, dt=1e-3, backend="numpy", tracker=tracker)
    e2 = EmulsionTimeCourse.from_file(path, progress=False)

    assert e1 == e2
    assert len(e1) == 2
    assert len(e1[0]) == 1  # found a single droplet
    assert path.stat().st_size > 0  # wrote some result

def test_length_scale_tracker(tmp_path):
    """test the length scale tracker"""
    path = tmp_path / "test_length_scale_tracker.json"

    grid = CartesianGrid([[0, 10 * np.pi]], 64, periodic=True)
    field = ScalarField.from_expression(grid, "sin(2 * x)")

    pde = DiffusionPDE()
    tracker = LengthScaleTracker(0.05, filename=path)
    pde.solve(field, t_range=0.1, backend="numpy", tracker=tracker)

    for ls in tracker.length_scales:
        assert ls == pytest.approx(np.pi, rel=1e-3)
    assert path.stat().st_size > 0  # wrote some result

max_stars_repo_name: tefavidal/py-droplets

max_issues_repo_name: tefavidal/py-droplets

max_forks_repo_name: tefavidal/py-droplets


text:
"""
Gillespie Stochastic Simulation Algorithm
----

This part of the package provides a simple implementation of GSSA.
This is designed for experimental purposes much more than for performance.
If you would like to use SSA for parameter inference, or for high
number of species, there are many superior implementations available.
"""

import multiprocessing
import numpy as np
from means.simulation.trajectory import Trajectory, TrajectoryCollection
from means.io.serialise import SerialisableObject
from means.util.moment_counters import generate_n_and_k_counters
from means.util.sympyhelpers import product, substitute_all
from means.core import Moment

class SSASimulation(SerialisableObject):
    """
        A class providing an implementation of the exact Gillespie Stochastic Simulation Algorithm [Gillespie77].

            >>> from means.examples import MODEL_P53
            >>> from means import StochasticProblem, SSASimulation
            >>> import numpy as np
            >>> PROBLEM = StochasticProblem(MODEL_P53)
            >>> RATES = [90, 0.002, 1.7, 1.1, 0.93, 0.96, 0.01]
            >>> INITIAL_CONDITIONS = [70, 30, 60]
            >>> TIME_RANGE = np.arange(0, 40, .1)
            >>> N_SSA = 10
            >>> ssas = SSASimulation(PROBLEM, N_SSA)
            >>> mean_trajectories = ssas.simulate_system(RATES, INITIAL_CONDITIONS, TIME_RANGE)

    .. [Gillespie77]Gillespie, Daniel T. "Exact stochastic simulation of coupled chemical reactions."\
         The journal of physical chemistry 81.25 (1977): 2340-2361.
    """
    def __init__(self, stochastic_problem, n_simulations, random_seed=None):
        """

        :param stochastic_problem:
        :param n_simulations:
        :param random_seed:
        """
        self.__random_seed = random_seed
        self.__problem = stochastic_problem
        self.__n_simulations = n_simulations

    def _validate_parameters(self, parameters, initial_conditions):

        if len(self.__problem.species) != len(initial_conditions):
            exception_str = "The number of initial conditions and the number of species are different. ({0} != {1})"
            raise Exception(exception_str.format(len(self.__problem.species), len(initial_conditions)))

        elif len(self.__problem.parameters) != len(parameters):
            exception_str = "The number of parameters and the number of constants are different. ({0} != {1})"
            raise Exception(exception_str.format(len(self.__problem.parameters), len(parameters)))

    def simulate_system(self, parameters, initial_conditions, timepoints,
                        max_moment_order=1, number_of_processes=1):
        """
        Perform Gillespie SSA simulations and returns trajectories for of each species.
        Each trajectory is interpolated at the given time points.
        By default, the average amounts of species for all simulations is returned.

        :param parameters: list of the initial values for the constants in the model.
                                  Must be in the same order as in the model
        :param initial_conditions: List of the initial values for the equations in the problem.
                        Must be in the same order as these equations occur.

        :param timepoints: A list of time points to simulate the system for

        :param number_of_processes: the number of parallel process to be run
        :param max_moment_order: the highest moment order to calculate the trajectories to.
                                 if set to zero, the individual trajectories will be returned, instead of
                                 the averaged moments.
        E.g. a value of one will return means, a values of two, means, variances and covariance and so on.

        :return: a list of :class:`~means.simulation.Trajectory` one per species in the problem,
            or a list of lists of trajectories (one per simulation) if `return_average == False`.
        :rtype: list[:class:`~means.simulation.Trajectory`]
        """
        max_moment_order = int(max_moment_order)
        assert(max_moment_order >= 0)

        n_simulations = self.__n_simulations
        self._validate_parameters(parameters, initial_conditions)
        t_max= max(timepoints)

        substitution_pairs = dict(zip(self.__problem.parameters, parameters))
        propensities = substitute_all(self.__problem.propensities, substitution_pairs)
        # lambdify the propensities for fast evaluation
        propensities_as_function = self.__problem.propensities_as_function
        def f(*species_parameters):
            return propensities_as_function(*(np.concatenate((species_parameters, parameters))))

        population_rates_as_function = f

        if not self.__random_seed:
            seed_for_processes = [None] * n_simulations
        else:
            seed_for_processes = [i for i in range(self.__random_seed, n_simulations + self.__random_seed)]

        if number_of_processes ==1:
            ssa_generator = _SSAGenerator(population_rates_as_function,
                                        self.__problem.change, self.__problem.species,
                                        initial_conditions, t_max, seed=self.__random_seed)

            results = map(ssa_generator.generate_single_simulation, seed_for_processes)

        else:
            p = multiprocessing.Pool(number_of_processes,
                    initializer=multiprocessing_pool_initialiser,
                    initargs=[population_rates_as_function, self.__problem.change,
                              self.__problem.species,
                              initial_conditions, t_max, self.__random_seed])

            results = p.map(multiprocessing_apply_ssa, seed_for_processes)

            p.close()
            p.join()

        resampled_results = [[traj.resample(timepoints, extrapolate=True) for traj in res] for res in results]
        for i in resampled_results:
            idx = len(i[0].values) - 1

        if max_moment_order == 0:
            # Return a list of TrajectoryCollection objects
            return map(TrajectoryCollection, resampled_results)

        moments = self._compute_moments(resampled_results, max_moment_order)
        return TrajectoryCollection(moments)

    def _compute_moments(self, all_trajectories, max_moment_order):

        mean_trajectories = [sum(trajs)/float(len(trajs)) for trajs in zip(*all_trajectories)]

        if max_moment_order == 1:
            return mean_trajectories
        n_counter, _ = generate_n_and_k_counters(max_moment_order - 1, self.__problem.species)

        out_trajects = mean_trajectories[:]
        for n in n_counter:
            if n.order == 0:
                continue
            out_trajects.append(self._compute_one_moment(all_trajectories, mean_trajectories, n))
        return out_trajects

    def _compute_one_moment(self, all_trajectories, mean_trajectories, moment):

        # the expectation of the product:
        #products_of_sps = [product(trajs) for trajs in all_trajectories]
        n_vec = moment.n_vector

        to_multipl = []

        for i, trajs in enumerate(zip(*all_trajectories)):
            mean_of_sp = mean_trajectories[i]
            order_of_sp = n_vec[i]
            xi_minus_ex = [(t - mean_of_sp) ** order_of_sp for t in trajs]
            for x in xi_minus_ex:
                x.set_description(moment)
            to_multipl.append(xi_minus_ex)

        to_sum = [product(xs) for xs in zip(*to_multipl)]

        return sum(to_sum)/ float(len(to_sum))

def multiprocessing_pool_initialiser(population_rates_as_function, change, species,
                                     initial_conditions, t_max, seed):
    global ssa_generator
    current = multiprocessing.current_process()
    #increment the random seed inside each process at creation, so the result should be reproducible
    if seed:
        seed += current._identity[0]
    else:
        seed = current._identity[0]
    ssa_generator = _SSAGenerator(population_rates_as_function, change, species, initial_conditions, t_max, seed)

def multiprocessing_apply_ssa(x):
    """
    Used in the SSASimulation class.
    Needs to be in global scope for multiprocessing module to pick it up
    """
    result = ssa_generator.generate_single_simulation(x)
    return result

class _SSAGenerator(object):
    def __init__(self, population_rates_as_function, change, species, initial_conditions, t_max, seed):
        """
        :param population_rates_as_function: function to evaluate propensities given the amount of species
        :param change: the change matrix (transpose of the stoichiometry matrix) as an numpy in array
        :param initial_conditions: the initial conditions of the system
        :param t_max: the time when the simulation should stop
        :param seed: an integer to initialise the random seed. If `None`, the random seed will be set
                automatically (e.g. from /dev/random) once for all.
        """
        self.__rng = np.random.RandomState(seed)
        self.__population_rates_as_function = population_rates_as_function
        self.__change = change
        self.__initial_conditions = initial_conditions
        self.__t_max = t_max
        self.__species = species

    def _gssa(self, initial_conditions, t_max):
        """
        This function is inspired from Yoav Ram's code available at:
        http://nbviewer.ipython.org/github/yoavram/ipython-notebooks/blob/master/GSSA.ipynb

        :param initial_conditions: the initial conditions of the system
        :param t_max:  the time when the simulation should stop
        :return:
        """
        # set the initial conditions and t0 = 0.
        species_over_time = [np.array(initial_conditions).astype("int16")]
        t = 0
        time_points = [t]
        while t < t_max and species_over_time[-1].sum() > 0:
            last = species_over_time[-1]
            e, dt = self._draw(last)
            t += dt
            species_over_time.append(last + self.__change[e,:])
            time_points.append(t)
        return time_points, np.array(species_over_time).T

    def _draw(self, population):
        population_rates = np.array(self.__population_rates_as_function(*population))
        sum_rates = population_rates.sum()
        time = self.__rng.exponential(1.0/sum_rates)
        event_distribution = population_rates.flatten()
        event_distribution /= event_distribution.sum()
        event = self.__rng.multinomial(1, event_distribution).argmax()
        return event, time

    def generate_single_simulation(self, x):
        """
        Generate a single SSA simulation
        :param x: an integer to reset the random seed. If None, the initial random number generator is used
        :return: a list of :class:`~means.simulation.Trajectory` one per species in the problem
        :rtype: list[:class:`~means.simulation.Trajectory`]
        """
        #reset random seed
        if x:
            self.__rng = np.random.RandomState(x)

        # perform one stochastic simulation
        time_points, species_over_time = self._gssa(self.__initial_conditions, self.__t_max)

        # build descriptors for first order raw moments aka expectations (e.g. [1, 0, 0], [0, 1, 0] and [0, 0, 1])
        descriptors = []
        for i, s in enumerate(self.__species):
            row = [0] * len(self.__species)
            row[i] = 1
            descriptors.append(Moment(row, s))

        # build trajectories
        trajectories = [Trajectory(time_points, spot, desc) for
                        spot, desc in zip(species_over_time, descriptors)]

        return trajectories

max_stars_repo_name: nicktimko/means

max_issues_repo_name: nicktimko/means

max_forks_repo_name: nicktimko/means


text:
import argparse
import math

import numpy as np

EXAMPLE_USAGE = """
Example Usage:
    python upper_bound --n_agents 4 --episode_len 1000 --n_apples 20
"""

def create_parser(parser_creator=None):
    parser_creator = parser_creator or argparse.ArgumentParser
    parser = parser_creator(
        formatter_class=argparse.RawDescriptionHelpFormatter,
        description="Determine maximum number of gathered apples",
        epilog=EXAMPLE_USAGE)

    required_named = parser.add_argument_group("required named arguments")
    required_named.add_argument(
        "--n_agents",
        type=int,
        required=True,
        help="Number of agents")
    required_named.add_argument(
        "--episode_len",
        type=int,
        required=True,
        help="Length of the episode")
    required_named.add_argument(
        "--n_apples",
        type=int,
        required=True,
        help="Initial number of spawned apples")
    required_named.add_argument(
        "--n_apple_frames",
        type=int,
        required=True,
        help="Number of frames after which apples are respawned")

    return parser

class UpperBoundSimulator:

    def __init__(self, n_apples, n_apple_frames, episode_length):
        self.apple_spawn_queue = np.zeros(n_apple_frames, dtype=np.uint8)
        self.max_n_apples_in_this_turn = n_apples
        self.n_apples = n_apples
        self.n_apples_init = n_apples
        self.episode_length = episode_length

    def n_apples_at_step(self, step_number):
        return round(
            ((math.cos(step_number * 2 * math.pi / self.episode_length) + 1) * 0.75 + 0.5) * self.n_apples_init / 2)

    def gather_and_enqueue_apples_to_respawn(self, step_number, apples_to_be_spawned):
        new_max_n_apples_in_this_turn = self.n_apples_at_step(step_number)
        if apples_to_be_spawned > 0 and self.max_n_apples_in_this_turn > new_max_n_apples_in_this_turn:
            # using min function to discount the remaining apples in the next steps
            difference = min(apples_to_be_spawned, self.max_n_apples_in_this_turn - new_max_n_apples_in_this_turn)
            self.max_n_apples_in_this_turn -= difference
            apples_to_be_spawned -= difference
        elif self.max_n_apples_in_this_turn < new_max_n_apples_in_this_turn:
            difference = new_max_n_apples_in_this_turn - self.max_n_apples_in_this_turn
            self.max_n_apples_in_this_turn += difference
            apples_to_be_spawned += difference

        self.queue_apples_to_be_spawned(apples_to_be_spawned)

    def spawn_apples(self):
        self.n_apples += self.apple_spawn_queue[0]
        self.apple_spawn_queue[0] = 0
        self.apple_spawn_queue = np.roll(self.apple_spawn_queue, -1)

    def queue_apples_to_be_spawned(self, apples: int):
        self.apple_spawn_queue[-1] += apples

if __name__ == "__main__":

    parser = create_parser()
    args = parser.parse_args()

    u_simulator = UpperBoundSimulator(args.n_apples, args.n_apple_frames, args.episode_len)

    consumed_total = 0

    for t in range(1, args.episode_len + 1):
        consumed = min(args.n_agents, u_simulator.n_apples)
        u_simulator.n_apples -= consumed
        consumed_total += consumed
        u_simulator.gather_and_enqueue_apples_to_respawn(t, consumed)
        u_simulator.spawn_apples()

    print(consumed_total, "apples consumed")


text:
import inviwopy as ivw
import ivwdataframe as df
import atomdata

import functools
import itertools
import math
import numpy

from pathlib import Path

def parseFile(file):
    file = Path(file)
    if file.suffix == ".xz":
        import lzma
        with lzma.open(file, mode='rt') as f:
            lines = f.readlines()
    elif file.suffix == ".bzip2" or file.suffix == ".bz2":
        import bz2
        with bz2.open(file, mode='rt') as f:
            lines = f.readlines()
    elif file.suffix == ".gz":
        import gzip
        with gzip.open(file, mode='rt') as f:
            lines = f.readlines()
    else:
        with open(file, mode='r') as f:
            lines = f.readlines()

    scale = float(lines[1].strip())
    a1 = numpy.array(list(map(float, lines[2].strip().split())))
    a2 = numpy.array(list(map(float, lines[3].strip().split())))
    a3 = numpy.array(list(map(float, lines[4].strip().split())))

    elem = lines[5].strip().split()
    nelem = list(map(int, lines[6].strip().split()))
    direct = not {lines[7].strip()[0]}.issubset({"K", "k", "C", "c"})
    ntot = functools.reduce(lambda x, y: x + y, nelem)
    pos = []
    for line in lines[8:8+ntot]:
        pos.append(list(map(float, line.strip().split())))
    pos = numpy.array(pos).astype(numpy.float32)

    factor = scale if scale >= 0.0 else numpy.power(
        -scale / numpy.dot(a3, numpy.cross(a2, a1)), 1/3)
    basis = numpy.array([a1, a2, a3]) * factor
    offset = -0.5 * (basis[0] + basis[1] + basis[2])
    if not direct:
        inv = numpy.linalg.inv(basis)
        pos = numpy.array([numpy.dot(inv, p) for p in pos])

    dims = list(map(int, lines[9+ntot].split()))
    voxels = functools.reduce(lambda x, y: x * y, dims)
    perrow = len(lines[10+ntot].strip().split())
    rows = math.ceil(voxels/perrow)

    chg = []
    for line in lines[10+ntot:10+ntot+rows]:
        for c in line.strip().split():
            chg.append(float(c))

    chgdata = numpy.array(chg).astype(numpy.float32)
    chgdata.shape = [dims[0], dims[1], dims[2]]

    volume = ivw.data.Volume(chgdata)
    volume.data = chgdata
    volume.dataMap.dataRange = ivw.glm.dvec2(chgdata.min(), chgdata.max())
    volume.dataMap.valueRange = volume.dataMap.dataRange

    volume.basis = ivw.glm.mat3(basis)
    volume.offset = ivw.glm.vec3(offset)

    elemtype = []
    for i, n in enumerate(nelem):
        for x in range(n):
            elemtype.append(elem[i])

    return (volume, pos, elem, nelem, elemtype)

def parseCubeFile(file):
    file = Path(file)
    if file.suffix == ".xz":
        import lzma
        with lzma.open(file, mode='rt') as f:
            lines = f.readlines()
    elif file.suffix == ".bzip2" or file.suffix == ".bz2":
        import bz2
        with bz2.open(file, mode='rt') as f:
            lines = f.readlines()
    elif file.suffix == ".gz":
        import gzip
        with gzip.open(file, mode='rt') as f:
            lines = f.readlines()
    else:
        with open(file, mode='r') as f:
            lines = f.readlines()

    nVal = 1
    extraLines = 0
    splitted = lines[2].strip().split()
    numAtoms = int(splitted[0])
    origin = numpy.array(list(map(float, splitted[1:4])))
    if numAtoms < 0:
        numAtoms = -numAtoms
        extraLines = 1
    if len(splitted) > 4:
        nVal = int(splitted[4])

    unitsInAngstrom = False
    bohrToAngstrom = 0.529177249
    splitted = lines[3].strip().split()
    sizeX = int(splitted[0])
    if sizeX < 0:
        sizeX = -sizeX
        unitsInAngstrom = True
    a1 = sizeX * numpy.array(list(map(float, splitted[1:4])))

    splitted = lines[4].strip().split()
    sizeY = abs(int(splitted[0]))
    a2 = sizeY * numpy.array(list(map(float, splitted[1:4])))

    splitted = lines[5].strip().split()
    sizeZ = abs(int(splitted[0]))
    a3 = sizeZ * numpy.array(list(map(float, splitted[1:4])))

    basis = numpy.array([a1, a2, a3]).astype(numpy.float32)
    if not unitsInAngstrom:
        basis = bohrToAngstrom * basis
        origin = bohrToAngstrom * origin
    offset = -0.5 * (basis[0] + basis[1] + basis[2])

    pos = []
    atomType = []
    for i in range(numAtoms):
        splitted = lines[6 + i].strip().split()
        atomID = int(splitted[0])
        atomPos = numpy.array(list(map(float, splitted[2:5])))
        if not unitsInAngstrom:
            atomPos = bohrToAngstrom * atomPos
        atomPos = atomPos - origin
        pos.append(atomPos)
        atomType.append(atomdata.atomicSymbol(atomID))

    chg = []
    for i in range(6 + numAtoms + extraLines, len(lines)):
        splitted = lines[i].strip().split()
        for val in splitted:
            chg.append(float(val))

    chgdata = numpy.array(chg).astype(numpy.float32)
    chosenIdx = 0
    chgdata = chgdata[chosenIdx::nVal]
    size = (sizeX, sizeY, sizeZ)
    chgdata = chgdata.reshape(size)
    chgdata = chgdata.flatten('F')
    chgdata = chgdata.reshape(size)

    volume = ivw.data.Volume(chgdata)
    volume.data = chgdata
    volume.dataMap.dataRange = ivw.glm.dvec2(chgdata.min(), chgdata.max())
    volume.dataMap.valueRange = volume.dataMap.dataRange

    volume.basis = ivw.glm.mat3(basis)
    volume.offset = ivw.glm.vec3(offset)

    return (volume, pos, atomType)

def createMesh(pos, elemtype, basis, offset, pm, margin):
    position = []
    color = []
    radius = []
    picking = []
    index = []

    pm.resize(len(elemtype))

    for i, p in enumerate(pos):
        c = atomdata.color(elemtype[i])
        r = atomdata.radius(elemtype[i])
        pi = pm.pickingId(i)

        def addVertex(vertexpos):
            position.append(vertexpos)
            color.append(c)
            radius.append(r)
            picking.append(pi)
            index.append(i)

        if margin > 0.0:
            for shift in itertools.product([-1, 0, 1], repeat=3):
                if all([x+s > -margin and x+s < 1 + margin for s, x in zip(shift, p)]):
                    addVertex([x+s for s, x in zip(shift, p)])
        else:
            addVertex(p)

    mesh = ivw.data.Mesh()
    mesh.basis = ivw.glm.mat3(basis)
    mesh.offset = ivw.glm.vec3(offset)

    mesh.addBuffer(ivw.data.BufferType.PositionAttrib, ivw.data.Buffer(
        numpy.array(position).astype(numpy.float32)))
    mesh.addBuffer(ivw.data.BufferType.ColorAttrib, ivw.data.Buffer(
        numpy.array(color).astype(numpy.float32)))
    mesh.addBuffer(ivw.data.BufferType.RadiiAttrib, ivw.data.Buffer(
        numpy.array(radius).astype(numpy.float32)))
    mesh.addBuffer(ivw.data.BufferType.PickingAttrib, ivw.data.Buffer(
        numpy.array(picking).astype(numpy.uint32)))
    mesh.addBuffer(ivw.data.BufferType.IndexAttrib, ivw.data.Buffer(
        numpy.array(index).astype(numpy.uint32)))
    return mesh

def createMeshForCube(pos, elemtype, basis, offset, pm):
    position = []
    color = []
    radius = []
    picking = []
    index = []

    pm.resize(len(elemtype))

    for i, p in enumerate(pos):
        c = atomdata.color(elemtype[i])
        r = atomdata.radius(elemtype[i])
        pi = pm.pickingId(i)

        def addVertex(vertexpos):
            position.append(vertexpos)
            color.append(c)
            radius.append(r)
            picking.append(pi)
            index.append(i)

        addVertex(p)

    mesh = ivw.data.Mesh()
    #mesh.basis = ivw.glm.mat3(basis)
    mesh.offset = ivw.glm.vec3(offset)

    mesh.addBuffer(ivw.data.BufferType.PositionAttrib, ivw.data.Buffer(
        numpy.array(position).astype(numpy.float32)))
    mesh.addBuffer(ivw.data.BufferType.ColorAttrib, ivw.data.Buffer(
        numpy.array(color).astype(numpy.float32)))
    mesh.addBuffer(ivw.data.BufferType.RadiiAttrib, ivw.data.Buffer(
        numpy.array(radius).astype(numpy.float32)))
    mesh.addBuffer(ivw.data.BufferType.PickingAttrib, ivw.data.Buffer(
        numpy.array(picking).astype(numpy.uint32)))
    mesh.addBuffer(ivw.data.BufferType.IndexAttrib, ivw.data.Buffer(
        numpy.array(index).astype(numpy.uint32)))
    return mesh

def createDataFrame(pos, elemtype, modelMat):
    dataframe = df.DataFrame()
    ct = dataframe.addCategoricalColumn("type")
    cx = dataframe.addFloatColumn("x")
    cy = dataframe.addFloatColumn("y")
    cz = dataframe.addFloatColumn("z")

    for et, p in zip(elemtype, pos):
        mp = modelMat * ivw.glm.vec4(p[0], p[1], p[2], 1.0)

        ct.add(et)
        cx.add(mp[0])
        cy.add(mp[1])
        cz.add(mp[2])

    dataframe.updateIndex()

    return dataframe

def createDataFrameForCube(pos, elemtype):
    dataframe = df.DataFrame()
    ct = dataframe.addCategoricalColumn("type")
    cx = dataframe.addFloatColumn("x")
    cy = dataframe.addFloatColumn("y")
    cz = dataframe.addFloatColumn("z")
    r = dataframe.addFloatColumn("r")

    for et, p in zip(elemtype, pos):
        ct.add(et)
        cx.add(p[0])
        cy.add(p[1])
        cz.add(p[2])
        r.add(atomdata.radius(et))

    dataframe.updateIndex()

    return dataframe

max_stars_repo_name: grassofsky/modules

max_issues_repo_name: grassofsky/modules

max_forks_repo_name: grassofsky/modules


text:
from sc2.position import Point2
from sc2.unit import Unit
from sc2.units import Units
from sc2.game_state import GameState

import logging

logger = logging.getLogger(__name__)

import math
from math import pow
import numpy as np
import warnings

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    from scipy.spatial.distance import pdist, cdist

from typing import Dict, Tuple, Iterable, Generator

class DistanceCalculation:
    def __init__(self):
        self.state: GameState = None
        # self._generated_frame = -100
        self._generated_frame2 = -100
        # Pdist condensed vector generated by scipy pdist, half the size of the cdist matrix as 1d array
        self._cached_pdist: np.ndarray = None
        self._cached_cdist: np.ndarray = None

    @property
    def _units_count(self) -> int:
        return len(self.all_units)

    @property
    def _pdist(self) -> np.ndarray:
        """ As property, so it will be recalculated each time it is called, or return from cache if it is called multiple times in teh same game_loop. """
        if self._generated_frame2 != self.state.game_loop:
            return self.calculate_distances()
        return self._cached_pdist

    @property
    def _cdist(self) -> np.ndarray:
        """ As property, so it will be recalculated each time it is called, or return from cache if it is called multiple times in teh same game_loop. """
        if self._generated_frame2 != self.state.game_loop:
            return self.calculate_distances()
        return self._cached_cdist

    def _calculate_distances_method1(self) -> np.ndarray:
        self._generated_frame2 = self.state.game_loop
        # Converts tuple [(1, 2), (3, 4)] to flat list like [1, 2, 3, 4]
        flat_positions = (coord for unit in self.all_units for coord in unit.position_tuple)
        # Converts to numpy array, then converts the flat array back to shape (n, 2): [[1, 2], [3, 4]]
        positions_array: np.ndarray = np.fromiter(flat_positions, dtype=np.float, count=2 * self._units_count).reshape(
            (self._units_count, 2)
        )
        assert len(positions_array) == self._units_count
        # See performance benchmarks
        self._cached_pdist = pdist(positions_array, "sqeuclidean")

        return self._cached_pdist

    def _calculate_distances_method2(self) -> np.ndarray:
        self._generated_frame2 = self.state.game_loop
        # Converts tuple [(1, 2), (3, 4)] to flat list like [1, 2, 3, 4]
        flat_positions = (coord for unit in self.all_units for coord in unit.position_tuple)
        # Converts to numpy array, then converts the flat array back to shape (n, 2): [[1, 2], [3, 4]]
        positions_array: np.ndarray = np.fromiter(flat_positions, dtype=np.float, count=2 * self._units_count).reshape(
            (self._units_count, 2)
        )
        assert len(positions_array) == self._units_count
        # See performance benchmarks
        self._cached_cdist = cdist(positions_array, positions_array, "sqeuclidean")

        return self._cached_cdist

    def _calculate_distances_method3(self) -> np.ndarray:
        """ Nearly same as above, but without asserts"""
        self._generated_frame2 = self.state.game_loop
        flat_positions = (coord for unit in self.all_units for coord in unit.position_tuple)
        positions_array: np.ndarray = np.fromiter(flat_positions, dtype=np.float, count=2 * self._units_count).reshape(
            (-1, 2)
        )
        # See performance benchmarks
        self._cached_cdist = cdist(positions_array, positions_array, "sqeuclidean")

        return self._cached_cdist

    # Helper functions

    def square_to_condensed(self, i, j) -> int:
        # Converts indices of a square matrix to condensed matrix
        # https://stackoverflow.com/a/36867493/10882657
        assert i != j, "No diagonal elements in condensed matrix! Diagonal elements are zero"
        if i < j:
            i, j = j, i
        return self._units_count * j - j * (j + 1) // 2 + i - 1 - j

    def convert_tuple_to_numpy_array(self, pos: Tuple[float, float]) -> np.ndarray:
        """ Converts a single position to a 2d numpy array with 1 row and 2 columns. """
        return np.fromiter(pos, dtype=float, count=2).reshape((1, 2))

    # Fast and simple calculation functions

    def distance_math_hypot(self, p1: Tuple[float, float], p2: Tuple[float, float]):
        return math.hypot(p1[0] - p2[0], p1[1] - p2[1])

    def distance_math_hypot_squared(self, p1: Tuple[float, float], p2: Tuple[float, float]):
        return pow(p1[0] - p2[0], 2) + pow(p1[1] - p2[1], 2)

    def _distance_squared_unit_to_unit_method0(self, unit1: Unit, unit2: Unit) -> float:
        return self.distance_math_hypot_squared(unit1.position_tuple, unit2.position_tuple)

    # Distance calculation using the pre-calculated matrix above

    def _distance_squared_unit_to_unit_method1(self, unit1: Unit, unit2: Unit) -> float:
        # If checked on units if they have the same tag, return distance 0 as these are not in the 1 dimensional pdist array - would result in an error otherwise
        if unit1.tag == unit2.tag:
            return 0
        # Calculate index, needs to be after pdist has been calculated and cached
        condensed_index = self.square_to_condensed(unit1.distance_calculation_index, unit2.distance_calculation_index)
        assert condensed_index < len(
            self._cached_pdist
        ), f"Condensed index is larger than amount of calculated distances: {condensed_index} < {len(self._cached_pdist)}, units that caused the assert error: {unit1} and {unit2}"
        distance = self._pdist[condensed_index]
        return distance

    def _distance_squared_unit_to_unit_method2(self, unit1: Unit, unit2: Unit) -> float:
        # Calculate index, needs to be after cdist has been calculated and cached
        return self._cdist[unit1.distance_calculation_index, unit2.distance_calculation_index]

    # Distance calculation using the fastest distance calculation functions

    def _distance_pos_to_pos(self, pos1: Tuple[float, float], pos2: Tuple[float, float]) -> float:
        return self.distance_math_hypot(pos1, pos2)

    def _distance_units_to_pos(self, units: Units, pos: Tuple[float, float]) -> Generator[float, None, None]:
        """ This function does not scale well, if len(units) > 100 it gets fairly slow """
        return (self.distance_math_hypot(u.position_tuple, pos) for u in units)

    def _distance_unit_to_points(
        self, unit: Unit, points: Iterable[Tuple[float, float]]
    ) -> Generator[float, None, None]:
        """ This function does not scale well, if len(points) > 100 it gets fairly slow """
        pos = unit.position_tuple
        return (self.distance_math_hypot(p, pos) for p in points)

    def _distances_override_functions(self, method: int = 0):
        """ Overrides the internal distance calculation functions at game start in bot_ai.py self._prepare_start() function
        method 0: Use python's math.hypot
        The following methods calculate the distances between all units once:
        method 1: Use scipy's pdist condensed matrix (1d array)
        method 2: Use scipy's cidst square matrix (2d array)
        method 3: Use scipy's cidst square matrix (2d array) without asserts (careful: very weird error messages, but maybe slightly faster) """
        assert 0 <= method <= 3, f"Selected method was: {method}"
        if method == 0:
            self._distance_squared_unit_to_unit = self._distance_squared_unit_to_unit_method0
        elif method == 1:
            self._distance_squared_unit_to_unit = self._distance_squared_unit_to_unit_method1
            self.calculate_distances = self._calculate_distances_method1
        elif method == 2:
            self._distance_squared_unit_to_unit = self._distance_squared_unit_to_unit_method2
            self.calculate_distances = self._calculate_distances_method2
        elif method == 3:
            self._distance_squared_unit_to_unit = self._distance_squared_unit_to_unit_method2
            self.calculate_distances = self._calculate_distances_method3


text:
# -*- coding: utf-8 -*-
"""huckel_solver.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1E3qBURf9CByVwHQlv5IrV4shAxh-WA_w
"""

import numpy as np
from collections import Counter

# Returns list of eigenvalues for any square matrix
def get_evals(A: np.ndarray):
    evals, evecs = np.linalg.eig(A)
    # Rounding evals to avoid numerical instability
    evals_rounded = [round(eval, 5) for eval in evals]
    return evals_rounded

# Returns n x n Huckel matrix for a linear n polyene, where n >= 2
def get_H_linear(n: int):
    assert(n >= 2)
    H = np.zeros((n, n))
    # Hardcoding the first and last rows
    H[0, 1] = -1
    H[n-1, n-2] = -1
    if n != 2:
      i = 1
      while i < n-1:
        H[i, i-1] = -1
        H[i, i+1] = -1
        i += 1
    return H

# Returns n x n Huckel matrix for a cyclic n polyene, where n >= 3
def get_H_cyclic(n : int):
    assert(n >= 3)
    H = np.zeros((n, n))
    # Hardcoding the first and last rows
    H[0, 1], H[0, n-1] = -1, -1
    H[n-1, 0], H[n-1, n-2] = -1, -1
    i = 1
    while i < n-1:
      H[i, i-1] = -1
      H[i, i+1] = -1
      i += 1
    return H

# Returns n x n Huckel matrix for a solid (the five Platonic solids and Buckminsterfullerene)
def get_H_solid(n: int):
    if n not in [4, 6, 8, 12, 20, 60]:
        raise ValueError('Value of n does not correspond to a valid Platonic solid.')

    # Tetrahedron
    if n == 4:
        H = np.array([[0,-1,-1,-1],
                      [-1,0,-1,-1],
                      [-1,-1,0,-1],
                      [-1,-1,-1,0]])

    # Octahedron
    if n == 6:
        H = np.array([[0,-1,-1,-1,-1,0],
                      [-1,0,-1,-1,0,-1],
                      [-1,-1,0,0,-1,-1],
                      [-1,-1,0,0,-1,-1],
                      [-1,0,-1,-1,0,-1],
                      [0,-1,-1,-1,-1,0]])

    # Cube
    if n == 8:
        H = np.array([[0,0,0,0,0,-1,-1,-1],
                      [0,0,0,0,-1,0,-1,-1],
                      [0,0,0,0,-1,-1,0,-1],
                      [0,0,0,0,-1,-1,-1,0],
                      [0,-1,-1,-1,0,0,0,0],
                      [-1,0,-1,-1,0,0,0,0],
                      [-1,-1,0,-1,0,0,0,0],
                      [-1,-1,-1,0,0,0,0,0]])

    # Icosahedron
    if n == 12:
        H = np.array([[0,-1,0,0,-1,0,-1,-1,0,0,-1,0],
                      [-1,0,0,0,-1,-1,0,-1,0,-1,0,0],
                      [0,0,0,-1,0,0,-1,0,-1,0,-1,-1],
                      [0,0,-1,0,0,0,-1,-1,0,-1,0,-1],
                      [-1,-1,0,0,0,-1,0,0,-1,0,-1,0],
                      [0,-1,0,0,-1,0,0,0,-1,-1,0,-1],
                      [-1,0,-1,-1,0,0,0,-1,0,0,-1,0],
                      [-1,-1,0,-1,0,0,-1,0,0,-1,0,0],
                      [0,0,-1,0,-1,-1,0,0,0,0,-1,-1],
                      [0,-1,0,-1,0,-1,0,-1,0,0,0,-1],
                      [-1,0,-1,0,-1,0,-1,0,-1,0,0,0],
                      [0,0,-1,-1,0,-1,0,0,-1,-1,0,0]])

    # Dodecahedron
    if n == 20:
        H = np.array([[0,0,-1,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,-1,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0],
                      [-1,0,0,0,0,0,0,0,-1,-1,0,0,0,0,0,0,0,0,0,0],
                      [-1,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0],
                      [-1,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0],
                      [0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,-1,0,0,0,0],
                      [0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0],
                      [0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1],
                      [0,0,-1,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,-1,0,0],
                      [0,0,-1,0,0,0,0,0,0,0,0,-1,0,0,0,0,-1,0,0,0],
                      [0,0,0,-1,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,-1],
                      [0,0,0,0,-1,0,0,0,0,-1,0,0,0,0,0,0,0,0,-1,0],
                      [0,0,0,-1,0,0,0,0,0,0,0,0,0,-1,-1,0,0,0,0,0],
                      [0,0,0,0,-1,0,0,0,0,0,0,0,-1,0,0,-1,0,0,0,0],
                      [0,0,0,0,0,-1,0,0,0,0,0,0,-1,0,0,0,0,0,0,-1],
                      [0,0,0,0,0,-1,0,0,0,0,0,0,0,-1,0,0,0,0,-1,0],
                      [0,0,0,0,0,0,-1,0,0,-1,0,0,0,0,0,0,0,-1,0,0],
                      [0,0,0,0,0,0,0,-1,-1,0,0,0,0,0,0,0,-1,0,0,0],
                      [0,0,0,0,0,0,-1,0,0,0,0,-1,0,0,0,-1,0,0,0,0],
                      [0,0,0,0,0,0,0,-1,0,0,-1,0,0,0,-1,0,0,0,0,0]])

    # Buckminsterfullerene
    if n == 60:
        H = np.array([[0,-1,0,0,-1,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [-1,0,-1,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [-1,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [-1,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,-1,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,-1,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,-1,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,-1,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0,-1],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0,-1],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,-1,0,-1,0],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,-1,0,-1],
                      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,-1,0,0,-1,0]])

    return H

def get_eval_degeneracies(evals: list):
    c = Counter(evals)
    return c

def test(n):
    print(get_eval_degeneracies(get_evals(get_H_linear(n))))
    print(get_eval_degeneracies(get_evals(get_H_cyclic(n))))
    for i in [4, 6, 8, 12, 20, 60]:
        print(get_eval_degeneracies(get_evals(get_H_solid(i))))

if __name__ == '__main__':
    mol_type = int(input())

    # 0 corresponds to linear
    if mol_type == 0:
        print(get_eval_degeneracies(get_evals(get_H_linear(int(input())))))

    # 1 corresponds to cyclic
    if mol_type == 1:
        print(get_eval_degeneracies(get_evals(get_H_cyclic(int(input())))))

    if mol_type == 2:
        print(get_eval_degeneracies(get_evals(get_H_solid(int(input())))))

max_stars_repo_name: DavidD-Chen/huckel-solver

max_issues_repo_name: DavidD-Chen/huckel-solver

max_forks_repo_name: DavidD-Chen/huckel-solver


text:
# -*- coding: utf-8 -*-
r""".. _sphara_basis_eeg:

Determination of the SPHARA basis functions for an EEG sensor setup
===================================================================

.. topic:: Section contents

   This tutorial introduces the steps necessary to determine a
   generalized spatial Fourier basis for an :term:`EEG` sensor setup
   using SpharaPy. The special properties of the different
   discretization approaches of the Laplace-Beltrami operator will be
   discussed.

Introduction
------------

A Fourier basis is a solution to Laplace's eigenvalue problem

.. math::

   L \vec{x} = \lambda \vec{x}\,,\qquad\qquad(1)

with the discrete :term:`Laplace-Beltrami operator` in matrix notation
:math:`L \in \mathbb{R}^{M \times N}`, the eigenvectors
:math:`\vec{x}` containing the harmonic functions and the eigenvalues
:math:`\lambda` the natural frequencies.

By solving a Laplace eigenvalue problem, it is also possible to
determine a basis for a spatial Fourier analysis. Often in practical
applications, a measured quantity to be subjected to a Fourier
analysis is only known at spatially discrete sampling points (the
sensor positions). An arbitrary arrangement of sample points on a
surface in three-dimensional space can be described by means of a
:term:`triangular mesh`. In the case of an :term:`EEG` system, the
sample positions (the vertices of the triangular mesh) are the
locations of the sensors arranged on the head surface. A SPHARA basis
is a solution of a Laplace eigenvalue problem for the given triangle
mesh, that can be obtained by discretizing a Laplace-Beltrami operator
for the mesh and solving the Laplace eigenvalue problem in equation
(1). The SpharaPy package provides three methods for the
discretization of the Laplace-Beltrami operator; unit weighting of the
edges and weighting with the inverse of the Euclidean distance of the
edges, and a FEM approach. For more detailed information please refer
to :ref:`introduction` and :ref:`eigensystems_of_lb_operators`.

"""

######################################################################
# At the beginning we import three modules of the SpharaPy package as
# well as several other packages and single functions of packages.

# Code source: Uwe Graichen
# License: BSD 3 clause

# import modules from spharapy package
import spharapy.trimesh as tm
import spharapy.spharabasis as sb
import spharapy.datasets as sd

# import additional modules used in this tutorial
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

######################################################################
# Specification of the spatial configuration of the EEG sensors
# -------------------------------------------------------------
#
# Import information about EEG sensor setup of the sample data set
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
#
# In this tutorial we will determine a SPHARA basis for a 256 channel
# :term:`EEG` system with equidistant layout. The data set is one of
# the example data sets contained in the SpharaPy toolbox, see
# :mod:`spharapy.datasets` and
# :func:`spharapy.datasets.load_eeg_256_channel_study`.

# loading the 256 channel EEG dataset from spharapy sample datasets
mesh_in = sd.load_eeg_256_channel_study()

######################################################################
# The dataset includes lists of vertices, triangles, and sensor
# labels, as well as :term:`EEG` data from previously performed
# experiment addressing the cortical activation related to
# somatosensory-evoked potentials (SEP).

print(mesh_in.keys())

######################################################################
# The triangulation of the :term:`EEG` sensor setup consists of 256
# vertices and 480 triangles.

vertlist = np.array(mesh_in['vertlist'])
trilist = np.array(mesh_in['trilist'])
print('vertices = ', vertlist.shape)
print('triangles = ', trilist.shape)

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

fig = plt.figure()
fig.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.02)
ax = fig.gca(projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('The triangulated EEG sensor setup')
ax.view_init(elev=20., azim=80.)
ax.set_aspect('auto')
ax.plot_trisurf(vertlist[:, 0], vertlist[:, 1], vertlist[:, 2],
                triangles=trilist, color='lightblue', edgecolor='black',
                linewidth=0.5, shade=True)
plt.show()

######################################################################
# Create a SpharaPy TriMesh instance
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
# In the next step we create an instance of the class
# :class:`spharapy.trimesh.TriMesh` from the list of vertices and
# triangles.

# create an instance of the TriMesh class
mesh_eeg = tm.TriMesh(trilist, vertlist)

######################################################################
# The class :class:`spharapy.trimesh.TriMesh` provides a
# number of methods to determine certain properties of the triangle
# mesh required to generate the SPHARA basis, listed below:

# print all implemented methods of the TriMesh class
print([func for func in dir(tm.TriMesh) if not func.startswith('__')])

######################################################################
# Determining SPHARA bases using different discretisation approaches
# ------------------------------------------------------------------
# Computing the basis functions
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
#
# In the final step of the tutorial we will calculate SPHARA bases for
# the given :term:`EEG` sensor setup. For this we create three
# instances of the class :class:`spharapy.spharabasis.SpharaBasis`. We
# use the three discretization approaches implemented in this class
# for the Laplace-Beltrami operator: unit weighting ('unit') and
# inverse Euclidean weigthing ('inv_euclidean') of the edges of the
# triangular mesh as well as the FEM discretization ('fem')

# 'unit' discretization
sphara_basis_unit = sb.SpharaBasis(mesh_eeg, 'unit')
basis_functions_unit, natural_frequencies_unit = sphara_basis_unit.basis()

# 'inv_euclidean' discretization
sphara_basis_ie = sb.SpharaBasis(mesh_eeg, 'inv_euclidean')
basis_functions_ie, natural_frequencies_ie = sphara_basis_ie.basis()

# 'fem' discretization
sphara_basis_fem = sb.SpharaBasis(mesh_eeg, 'fem')
basis_functions_fem, natural_frequencies_fem = sphara_basis_fem.basis()

######################################################################
# Visualization the basis functions
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
# The first 15 spatially low-frequency SPHARA basis functions are
# shown below, starting with DC at the top left.
#
# SPHARA basis using the discretization approache 'unit'
# """"""""""""""""""""""""""""""""""""""""""""""""""""""

# sphinx_gallery_thumbnail_number = 2
figsb1, axes1 = plt.subplots(nrows=5, ncols=3, figsize=(8, 12),
                             subplot_kw={'projection': '3d'})
for i in range(np.size(axes1)):
    colors = np.mean(basis_functions_unit[trilist, i + 0], axis=1)
    ax = axes1.flat[i]
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.view_init(elev=60., azim=80.)
    ax.set_aspect('auto')
    trisurfplot = ax.plot_trisurf(vertlist[:, 0], vertlist[:, 1],
                                  vertlist[:, 2], triangles=trilist,
                                  cmap=plt.cm.bwr,
                                  edgecolor='white', linewidth=0.)
    trisurfplot.set_array(colors)
    trisurfplot.set_clim(-0.15, 0.15)

cbar = figsb1.colorbar(trisurfplot, ax=axes1.ravel().tolist(), shrink=0.85,
                       orientation='horizontal', fraction=0.05, pad=0.05,
                       anchor=(0.5, -4.5))

plt.subplots_adjust(left=0.0, right=1.0, bottom=0.08, top=1.0)
plt.show()

######################################################################
# SPHARA basis using the discretization approache 'inv_euclidean'
# """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

figsb1, axes1 = plt.subplots(nrows=5, ncols=3, figsize=(8, 12),
                             subplot_kw={'projection': '3d'})
for i in range(np.size(axes1)):
    colors = np.mean(basis_functions_ie[trilist, i + 0], axis=1)
    ax = axes1.flat[i]
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.view_init(elev=60., azim=80.)
    ax.set_aspect('auto')
    trisurfplot = ax.plot_trisurf(vertlist[:, 0], vertlist[:, 1],
                                  vertlist[:, 2], triangles=trilist,
                                  cmap=plt.cm.bwr,
                                  edgecolor='white', linewidth=0.)
    trisurfplot.set_array(colors)
    trisurfplot.set_clim(-0.18, 0.18)

cbar = figsb1.colorbar(trisurfplot, ax=axes1.ravel().tolist(), shrink=0.85,
                       orientation='horizontal', fraction=0.05, pad=0.05,
                       anchor=(0.5, -4.5))

plt.subplots_adjust(left=0.0, right=1.0, bottom=0.08, top=1.0)
plt.show()

######################################################################
# SPHARA basis using the discretization approache 'fem'
# """""""""""""""""""""""""""""""""""""""""""""""""""""

figsb1, axes1 = plt.subplots(nrows=5, ncols=3, figsize=(8, 12),
                             subplot_kw={'projection': '3d'})
for i in range(np.size(axes1)):
    colors = np.mean(basis_functions_fem[trilist, i + 0], axis=1)
    ax = axes1.flat[i]
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.view_init(elev=60., azim=80.)
    ax.set_aspect('auto')
    trisurfplot = ax.plot_trisurf(vertlist[:, 0], vertlist[:, 1],
                                  vertlist[:, 2], triangles=trilist,
                                  cmap=plt.cm.bwr,
                                  edgecolor='white', linewidth=0.)
    trisurfplot.set_array(colors)
    trisurfplot.set_clim(-0.01, 0.01)

cbar = figsb1.colorbar(trisurfplot, ax=axes1.ravel().tolist(), shrink=0.85,
                       orientation='horizontal', fraction=0.05, pad=0.05,
                       anchor=(0.5, -4.5))

plt.subplots_adjust(left=0.0, right=1.0, bottom=0.08, top=1.0)
plt.show()


text:
from __future__ import division

import pytest

import numpy as np

from cudf.dataframe import Series

from . import utils

def test_series_ceil():
    arr = np.random.random(100) * 100
    sr = Series(arr)
    np.testing.assert_equal(sr.ceil().to_array(), np.ceil(arr))

def test_series_floor():
    arr = np.random.random(100) * 100
    sr = Series(arr)
    np.testing.assert_equal(sr.floor().to_array(), np.floor(arr))

@pytest.mark.parametrize('nelem', [1, 7, 8, 9, 32, 64, 128])
def test_validity_ceil(nelem):
    # Data
    data = np.random.random(nelem) * 100
    mask = utils.random_bitmask(nelem)
    bitmask = utils.expand_bits_to_bytes(mask)[:nelem]
    null_count = utils.count_zero(bitmask)
    sr = Series.from_masked_array(data, mask, null_count)

    # Result
    res = sr.ceil()

    na_value = -100000
    got = res.fillna(na_value).to_array()
    res_mask = np.asarray(bitmask, dtype=np.bool_)[:data.size]

    expect = np.ceil(data)
    expect[~res_mask] = na_value

    # Check
    print('expect')
    print(expect)
    print('got')
    print(got)

    np.testing.assert_array_equal(expect, got)