text:
import os
import sys
import pyhdfs
from pandas.core.algorithms import quantile
from pyspark.sql.types import StructType, StructField, StringType, IntegerType, FloatType, DoubleType
import pandas as pd
import numpy as np
import PinkMST
from Point import Point

from pyspark.sql.session import SparkSession
from pyspark import SparkConf
from pyspark.context import SparkContext
sc = SparkContext.getOrCreate(SparkConf().setMaster("local[*]").set("spark.cores.max","4"))

filePath = "test/data.txt"
hadoop_master ='hdfs://127.0.0.1:9000'

class DataSpliter():
    def __init__(self, inputFileName=filePath, numSplits=None, outputDir=None):

        self.inputFileName = inputFileName
        self.numSplits = numSplits
        self.outputDir = outputDir
        self.sc = sc
        self.data = None

    def getSparkConf(self):
        return self.sc

    def displayValue(points, size):
        print("=" * 100)
        num = 0
        for p in points:
            print("====:" + p)
            p += 1
            if p > size:
                break
        print("=" * 100)

    def readFile(self):
        data = self.sc.sequenceFile(self.outputDir).map(lambda x:x[1]).collect()
        PinkMST.displayValue(data,None,len(data))

    def createPartitionFiles(self,fileLoc,numPartitions):
        self.delete(fileLoc) #删除原有文件
        idSubgraphs = list(range(int(numPartitions)))
        self.sc.parallelize(idSubgraphs, numPartitions).saveAsTextFile("hdfs://127.0.0.1:9000" + fileLoc)

    def writeSequenceFile(self):
        self.sc.addPyFile('Point.py')
        data = self.loadData()
        numPoints = len(data)
        print("data number: ", numPoints)
        pointId = 0
        points = []
        for i in data:
            points.append(str(pointId) + ":" + str(i))
            pointId += 1

        self.displayValue(points, 5)
        self.saveAsSequenceFile(points)

    def displayValue(self,points, size):
        print("=" * 80)
        for i in range(size):
            print(points[i])
        print("=" * 80)

    def saveAsSequenceFile(self, points):
        self.delete(self.outputDir)
        pointsToWrite = self.sc.parallelize(points,self.numSplits)
        pointsPairToWrite = pointsToWrite.map(lambda x : x)
        pointsPairToWrite.saveAsTextFile(path=hadoop_master+self.outputDir)

    def loadData(self):
        data = self.sc.textFile(self.inputFileName)
        data = data.map(lambda line:line.split(" ")).map(lambda l:[float(i) for i in l])
        schema = StructType([
                StructField("f1", FloatType(), True),
                StructField("f2", FloatType(), True),
                StructField("f3", FloatType(), True),
                StructField("f4", FloatType(), True),
                StructField("f5", FloatType(), True),
#                 StructField("f6", FloatType(), True),
#                 StructField("f7", FloatType(), True),
#                 StructField("f8", FloatType(), True),
#                 StructField("f9", FloatType(), True),
#                 StructField("f10", FloatType(), True)
        ])
        data = SparkSession(self.sc).createDataFrame(data,schema)
        data = data.rdd.map(lambda line: [float(i) for i in line]).take(10000)
        return data

    def delete(self,file):
        fs = pyhdfs.HdfsClient(hosts='127.0.0.1,9000')
        if(fs.exists(file)):
            fs.delete(file,recursive=True)
        else:
            pass

    def writeSmallFiles(self):
        points = [Point(0, [0.0, 0.0]), Point(1, [1.0, 0.0]), Point(2, [0.0, 2.0]), Point(3, [3.0, 2.0])]
        self.saveAsSequenceFile_1(points, self.outputDir)

    @staticmethod
    def openSequenceFile(filename):
        return sc.sequenceFile(filename).collect()


text:
# Copyright 2021 Google LLC
#
# 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
#
#     https://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 os
import bem_2d
import time
import gmsh_parser

import numpy as np
import numpy as nnp
import pandas as pd
import scipy.special
import math
import matplotlib.pylab as pylab
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from matplotlib.collections import CircleCollection
import bem_3d
import backend
import matplotlib
from backend import np as np
import sys

from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from mpl_toolkits.mplot3d.art3d import Line3DCollection

def draw_surface(points, triangles, n, scalar, sample_points, frame, freq):
    fig = plt.figure()
    ax = Axes3D(fig=fig)
    min = nnp.min(scalar)
    max = nnp.max(scalar)
    scaled = nnp.minimum(nnp.maximum(0., (scalar - min) / (max - min)), 1.)
    colors = matplotlib.cm.jet(scaled)
    lighting = (nnp.abs(n[:, 1]))[..., None] * nnp.array(
        [1, 1, 1, 1])[None, :]  # n.l ish (directional light)
    lighting[:, 3] = 1  # set alpha to 1
    colors *= lighting

    raw = []
    lines = []
    for i in range(triangles.shape[0]):
        ps = [
            points[triangles[i, 0]], points[triangles[i, 1]],
            points[triangles[i, 2]]
        ]
        raw.append(nnp.stack(ps))
        c = .3333 * (ps[0] + ps[1] + ps[2])
    for sp in sample_points:
        for k in range(3):
            off = nnp.zeros(3)
            off[k] = .05
            lines.append((sp - off, sp + off))
    ax.view_init(elev=10., azim=80.)
    polys = Poly3DCollection(raw)
    edges = Line3DCollection(lines)
    edges.set_edgecolor((.5, 1, 1, 1.0))
    colors = nnp.maximum(0., nnp.minimum(colors, 1.))
    polys.set_facecolor(colors)
    polys.set_edgecolor((1, 1, 1, .3))
    ax.add_collection(polys)
    ax.add_collection(edges)
    ax.set_xlim(-.05, .05)
    ax.set_ylim(-.05, .05)
    ax.set_zlim(-.05, .05)
    if frame is not None:
        pylab.title('Freq %f Hz' % freq)
        pylab.savefig("frame.%d.png" % frame)

    else:
        pylab.show()

def run(input_mesh_filename):
    method = bem_3d.DUAL_SURFACE
    stuff = gmsh_parser.GmshParser(input_mesh_filename)
    points = stuff.nodes
    edges = stuff.triangles
    masks = stuff.triangle_physical_masks
    points /= 1000.  # convert from mm to meters
    emit = [x for x in masks.keys() if x.find("emit") != -1][0]
    mask = masks[emit]

    angles = np.linspace(-math.pi / 2, math.pi / 2, 65, dtype=np.float32)
    #print(angles)
    half = angles.shape[0] // 2
    num = angles.shape[0] // 2 // 4
    sub_angles = np.arange(half, angles.shape[0], num)

    mic_positions = np.array([
        np.sin(angles),
        np.cos(angles),
        np.zeros(angles.shape[0]),
    ])
    mic_positions_sub = nnp.array([
        np.sin(angles[sub_angles]),
        np.cos(angles[sub_angles]),
        np.zeros(sub_angles.shape[0])
    ])
    mic_positions_sub = np.transpose(mic_positions_sub)

    elems_x1 = points[edges[:, 0]]  # first point in element
    elems_x2 = points[edges[:, 1]]  # second point in element
    elems_x3 = points[edges[:, 2]]  # second point in element
    elem_centers = (.3333333333 * (elems_x1 + elems_x2 + elems_x3)).astype(
        np.float32)
    n, area = bem_3d.geom_v(elems_x1, elems_x2, elems_x3)
    N_elements = elems_x1.shape[0]
    if True:
        scalar = mask * 1.
        draw_surface(points, stuff.triangles, n, scalar, mic_positions_sub,
                     None, 0.)

    bins = np.logspace(np.log10(800), np.log10(20000), 128)

    solve_type = bem_3d.EXTERIOR
    datas = []
    for idx, frequency in enumerate(bins):
        with bem_3d.Timer("Whole solve"):
            print("method %s %d/%d freq %f" %
                  (method, idx, bins.shape[0], frequency))
            density = 1.29  # kg/m^3
            omega = 2 * math.pi * frequency  # Hz
            c = 344  # m/s
            k = bem_3d.WaveNumber(frequency, c)

            phi_w = np.zeros(N_elements)
            v_w = np.ones(N_elements)
            phi = np.zeros(N_elements, np.complex64)

            norm = n[:, 1]
            v = np.where(backend.index[mask],
                         norm * np.array(1. / (1.j * omega), np.complex64),
                         np.zeros(N_elements, np.complex64))
            phi, v = bem_3d.SolveHelmholtz_v(k,
                                             elems_x1,
                                             elems_x2,
                                             elems_x3,
                                             phi_w,
                                             phi,
                                             v_w,
                                             v,
                                             type=solve_type)
            # draw_surface(points, stuff.triangles, n, nnp.array(np.abs(phi)[:,0]), mic_positions_sub, frame=idx, freq=frequency)
            with bem_3d.Timer("Evaluate microphones"):
                vals = (bem_3d.EvaluatePosition_v(
                    k, elems_x1, elems_x2, elems_x3, phi, v,
                    np.transpose(mic_positions), solve_type) * omega * density)
                datas.append(vals)
    data = np.stack(datas, axis=1)

    stuff = []
    stuffnorm = []
    dbnorm = 20 * np.log10(np.abs(data / data[sub_angles[0], :][None, :]))
    db = 20 * np.log10(np.abs(data) / np.max(np.abs(data)))
    labels = []
    for idx in sub_angles:
        stuff.extend([bins, db[idx, :], '-'])
        stuffnorm.extend([bins, dbnorm[idx, :], '-'])
        labels.append(str(float(angles[idx] / math.pi * 180)))

    pylab.figure(figsize=(8, 16))
    pylab.subplot(3, 1, 1)
    pylab.set_cmap('jet')
    pylab.pcolor(db)
    pylab.colorbar()
    pylab.xlabel('frequency index')
    pylab.ylabel('angle index')
    pylab.subplot(3, 1, 2)
    pylab.plot(*stuff)
    pylab.title('Frequency Response')
    pylab.semilogx()
    pylab.xlabel('Freq Hz')
    pylab.ylabel('dB')
    pylab.legend(labels)

    pylab.subplot(3, 1, 3)

    pylab.plot(*stuffnorm)
    pylab.title('Frequency Response normalized to on-axis')
    pylab.semilogx()
    pylab.legend(labels)
    pylab.xlabel('Freq Hz')
    pylab.ylabel('dB')
    pylab.ylim([-30, 10])
    pylab.savefig('test.png')

# Check arguments
if len(sys.argv) != 2:
    print(f"Usage: {sys.argv[0]} <mesh file>")
    print(" Runs directivity simulation across frequency spectrum for mesh.")
    sys.exit(1)

# Run example
filename = sys.argv[1]
run(filename)

ext: py

max_stars_repo_name: aselle/wavextrema

max_issues_repo_name: aselle/wavextrema

max_forks_repo_name: aselle/wavextrema

reason: import numpy,import scipy


text:
import math
from copy import deepcopy
from functools import lru_cache
from typing import Iterable, List, Tuple

import numpy as np
import pyrr

from payton.math.matrix import IDENTITY_MATRIX, Matrix
from payton.math.vector import Vector3D

def to_4(v: Iterable, fill: float = 1.0) -> Vector3D:
    res = list(v)
    for _ in range(4 - len(res)):
        res.append(fill)
    return res

def sub_vector(v1: Vector3D, v2: Vector3D) -> Vector3D:
    """Substract vector V1 from vector V2 and return resulting Vector.

    Keyword arguments:
    v1 -- Vector 1
    v2 -- Vector 2
    """
    return [v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2]]

def cross_product(v1: Vector3D, v2: Vector3D) -> Vector3D:
    """Return the cross product vector of V1 and V2.

    Keyword arguments:
    v1 -- Vector 1
    v2 -- Vector 2
    """
    return [
        v1[1] * v2[2] - v1[2] * v2[1],
        v1[2] * v2[0] - v1[0] * v2[2],
        v1[0] * v2[1] - v1[1] * v2[0],
    ]

def vector_norm(v: Vector3D) -> float:
    """Calculate and return the Norm of a given Vector."""
    return math.sqrt((v[0] * v[0]) + (v[1] * v[1]) + (v[2] * v[2]))

def normalize_vector(v: Vector3D) -> Vector3D:
    """
    Normalize given vector and output normalized vector.

    Keyword Arguments:
    v -- Vector to normalize
    """
    vn = vector_norm(v)
    if vn == 0.0:
        return v
    invlen = 1.0 / vn
    return [v[0] * invlen, v[1] * invlen, v[2] * invlen]

def plane_normal(v1: Vector3D, v2: Vector3D, v3: Vector3D) -> Vector3D:
    """
    Calculate the normal of the plane defined by v1, v2, v3 Vectors in space.

    Keyword arguments:
    v1 -- V1 Vector of the plane (Corner of the triangle in the given plane)
    v2 -- V2 Vector of the plane (Corner of the triangle in the given plane)
    v3 -- V3 Vector of the plane (Corner of the triangle in the given plane)
    """
    return normalize_vector(cross_product(sub_vector(v2, v1), sub_vector(v3, v1)))

def vector_transform(v: Vector3D, matrix: Matrix) -> Vector3D:
    """
    Transform the given Vector into given Matrix coordinates.

    Keyword arguments:
    v -- Vector to transform
    matrix -- Matrix for the transformation
    """
    mx = matrix[0]
    my = matrix[1]
    mz = matrix[2]
    mw = matrix[3]

    rx = v[0] * mx[0] + v[1] * my[0] + v[2] * mz[0]
    ry = v[0] * mx[1] + v[1] * my[1] + v[2] * mz[1]
    rz = v[0] * mx[2] + v[1] * my[2] + v[2] * mz[2]

    return [rx + mw[0], ry + mw[1], rz + mw[2]]

def invert_vector(v: Vector3D) -> Vector3D:
    """
    Invert - Negate the given vector.

    Keyword arguments:
    v -- Vector to invert
    """
    return [-v[0], -v[1], -v[2]]

def min_max(vlist: List[Vector3D]) -> Tuple[Vector3D, Vector3D]:
    """
    Calculate the smallest and biggest components of the given vector list.

    Returns the result as Tuple of Vectors as minVector - maxVector

    Example:
        vlist = [
            Vector(-2, 3, 4),
            Vector(-3, -5, 9)
        ]

        min_max(vlist)

    would return:

        ([-3, -5, 4], [-2, 3, 9])

    Keyword arguments:
    vlist -- List of vectors
    """
    x_list = [item[0] for item in vlist]
    y_list = [item[1] for item in vlist]
    z_list = [item[2] for item in vlist]
    return (
        [min(x_list), min(y_list), min(z_list)],
        [max(x_list), max(y_list), max(z_list)],
    )

def scale_vector(v: Vector3D, factor: float) -> Vector3D:
    """
    Scale (multiply each component) of the given Vector and return as a result.

    Keyword arguments:
    v -- Vector to scale
    factor -- factor to multiply
    """
    return [v[0] * factor, v[1] * factor, v[2] * factor]

def add_vectors(v1: Vector3D, v2: Vector3D) -> Vector3D:
    """
    Add two vectors.

    Keyword arguments:
    v1 -- First vector
    v2 -- Second vector
    """
    return [v1[0] + v2[0], v1[1] + v2[1], v1[2] + v2[2]]

def dot_product(v1: Vector3D, v2: Vector3D) -> float:
    """
    Return vector dot product of v1 and v2.

    Keyword arguments:
    v1 -- First vector
    v2 -- Second vector
    """
    return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]

def vector_angle(v1: Vector3D, v2: Vector3D) -> float:
    """
    Calculate the angle between two given vectors.

    Keyword arguments:
    v1 -- First vector
    v2 -- Second vector
    """
    v1_n = normalize_vector(v1)
    v2_n = normalize_vector(v2)
    return math.acos(dot_product(v1_n, v2_n) / (vector_norm(v1_n) * vector_norm(v2_n)))

def mid_point(v1: Vector3D, v2: Vector3D) -> Vector3D:
    """
    Calculate the middle point between v1 and v2 vectors.

    Keyword arguments:
    v1 -- First vector
    v2 -- Second vector
    """
    diff = sub_vector(v2, v1)
    leng = vector_norm(diff)
    diff = normalize_vector(diff)
    diff = scale_vector(diff, leng / 2.0)
    return add_vectors(diff, v1)

def distance(v1: Vector3D, v2: Vector3D) -> float:
    """
    Calculate the distance between vectors v1 and v2.

    Keyword arguments:
    v1 -- First vector
    v2 -- Second vector
    """
    diff = sub_vector(v2, v1)
    return vector_norm(diff)

def rotate_around_z(v: Vector3D, angle: float) -> Vector3D:
    """Rotate the given vector around Z axis by given degrees.

    Keyword arguments:
    v -- Vector to rotate
    angle -- Angle of rotation in Radians
    """
    s = math.sin(angle)
    c = math.cos(angle)
    return [v[0] * c + v[1] * s, v[1] * c - v[0] * s, v[2]]

def create_rotation_matrix_raw(axis: Vector3D, angle: float) -> Matrix:
    """
    Create a rotation matrix as Numpy array.

    If you transform a Vector by the rotation Matrix, you will be rotating
    the given Vector according to the rotation defined by the matrix.

    If you multiply a matrix by the rotation matrix, resulting matrix will be rotated
    by the given rotation matrix.

    Keyword arguments:
    axis -- Rotation axis
    angle -- Angle of rotation in Radians
    """
    sin = math.sin(angle)
    cos = math.cos(angle)
    m_cos = 1 - cos
    axis = normalize_vector(axis)
    return [
        [
            (m_cos * axis[0] * axis[0]) + cos,
            (m_cos * axis[0] * axis[1]) - (axis[2] * sin),
            (m_cos * axis[2] * axis[0]) + (axis[1] * sin),
            0,
        ],
        [
            (m_cos * axis[0] * axis[1]) + (axis[2] * sin),
            (m_cos * axis[1] * axis[1]) + cos,
            (m_cos * axis[1] * axis[2]) - (axis[0] * sin),
            0,
        ],
        [
            (m_cos * axis[2] * axis[0]) - (axis[1] * sin),
            (m_cos * axis[1] * axis[2]) + (axis[0] * sin),
            (m_cos * axis[2] * axis[2]) + cos,
            0,
        ],
        [0.0, 0.0, 0.0, 1.0],
    ]

def create_rotation_matrix(axis: Vector3D, angle: float) -> np.ndarray:
    """
    Create rotation matrix as a numpy array.

    Numpy array is a well packed proper memory array that fits well into
    graphics card. Therefore, we are using Numpy Arrays to push stuff to
    the graphics card.

    If you transform a Vector by the rotation Matrix, you will be rotating
    the given Vector according to the rotation defined by the matrix.

    If you multiply a matrix by the rotation matrix, resulting matrix will be rotated
    by the given rotation matrix.

    Keyword arguments:
    axis -- Rotation axis
    angle -- Angle of rotation in Radians.
    """
    result = create_rotation_matrix_raw(axis, angle)
    return np.array(result, dtype=np.float32)

def scale_matrix(x: float, y: float, z: float) -> np.ndarray:
    """
    Create a scale matrix.

    If you multiply an objects matrix with a scale matrix, it will scale / resize
    the coordinates relative to the object. As a result, the object will look resized.

    Keyword arguments:
    x -- X Scale Factor
    y -- Y Scale Factor
    z -- Z Scale Factor
    """
    result = deepcopy(IDENTITY_MATRIX)
    result = [[x, 0, 0, 0], [0, y, 0, 0], [0, 0, z, 0], [0, 0, 0, 1]]
    return np.array(result, dtype=np.float32)

@lru_cache(maxsize=512)
def ortho(left: float, right: float, bottom: float, top: float) -> np.ndarray:
    """
    Create an orthographic projection matrix as a Numpy Array.

    Keyword arguments:
    left -- Viewport left side (-x) coordinate
    right -- Viewport right side (+x) coordinate
    bottom -- Viewport bottom side (+y) coordinate
    top -- Viewport top side (-y) coordinate
    """
    result = [
        [2 / (right - left), 0, 0, 0],
        [0, 2 / (top - bottom), 0, 0],
        [0, 0, -1, 0],
        [-(right + left) / (right - left), -(top + bottom) / (top - bottom), 0, 1],
    ]
    return np.array(result, dtype=np.float32)

def cubemap_projection_matrices(from_point: Vector3D, far_plane: float) -> List[np.ndarray]:
    """
    Create the required Cubemap projection matrices.

    This method is suitable for generating a Shadow Map.

    Simply speaking, this method generates 6 different camera matrices from the center of
    an imaginary cube and covers all surfaces without conflicting.

    Keyword arguments;
    from_point -- Imaginary camera location
    far_plane -- How far the camera is capable of seeing. (Effects performance!)
    """

    def a2np(a: List[float]) -> np.ndarray:
        return np.array(a, dtype=np.float32)

    shadow_proj = pyrr.matrix44.create_perspective_projection(90.0, 1.0, 0.01, far_plane, np.float32)
    lightpos = np.array(list(from_point), dtype=np.float32)[:3]

    nx = pyrr.matrix44.create_look_at(
        lightpos,
        np.array(
            lightpos + a2np([-1.0, 0, 0]),
            dtype=np.float32,
        ),
        a2np([0, -1.0, 0]),
        dtype=np.float32,
    )
    px = pyrr.matrix44.create_look_at(
        lightpos,
        np.array(
            lightpos + a2np([1, 0, 0]),
            dtype=np.float32,
        ),
        a2np([0, -1.0, 0]),
        dtype=np.float32,
    )
    ny = pyrr.matrix44.create_look_at(
        lightpos,
        np.array(
            lightpos + a2np([0, -1, 0]),
            dtype=np.float32,
        ),
        a2np([0, 0, -1.0]),
        dtype=np.float32,
    )
    py = pyrr.matrix44.create_look_at(
        lightpos,
        np.array(
            lightpos + a2np([0, 1, 0]),
            dtype=np.float32,
        ),
        a2np([0, 0, 1.0]),
        dtype=np.float32,
    )
    pz = pyrr.matrix44.create_look_at(
        lightpos,
        np.array(
            lightpos + a2np([0, 0, 1]),
            dtype=np.float32,
        ),
        a2np([0, -1.0, 0]),
        dtype=np.float32,
    )
    nz = pyrr.matrix44.create_look_at(
        lightpos,
        np.array(
            lightpos + a2np([0, 0, -1]),
            dtype=np.float32,
        ),
        a2np([0, -1.0, 0]),
        dtype=np.float32,
    )

    return [
        px.dot(shadow_proj),
        nx.dot(shadow_proj),
        py.dot(shadow_proj),
        ny.dot(shadow_proj),
        pz.dot(shadow_proj),
        nz.dot(shadow_proj),
    ]

max_stars_repo_name: sinanislekdemir/payton

max_issues_repo_name: sinanislekdemir/payton

max_forks_repo_name: sinanislekdemir/payton

reason: import numpy


text:
import concurrent.futures
import numpy as np
import argparse
import os

from copy import deepcopy
from time import time
from pprint import pprint
from evolution.operators import *
from to.probabilistic_model import ProbabilisticModel
from to.mixture_model import MixtureModel
from evolution.chromosome import *
from utils.double_pole_physics import PoledCart
from utils.neural_network import Net

def evolutionary_algorithm(problem, dims, psize=100, gen=100, src_models=[],
                           stop_condition=True, create_model=True, multi_proc=False,
                           workers=1):

  fitness_hist = np.zeros((gen, psize))
  fitness_time = np.zeros((gen))

  init_func = lambda n: np.round(np.random.rand(n))
  pop = get_pop_init(psize, dims, init_func)
  start = time()

  if multi_proc:
    with concurrent.futures.ProcessPoolExecutor(max_workers=workers) as executor:
        processes = []
        for j in range(psize):

            proc = executor.submit(pop[j].fitness_calc, problem)
            processes.append(proc)

        fitnesses = [proc.result() for proc in processes]

        for j in range(psize):
            pop[j].fitness = fitnesses[j]
  else:
    fitnesses = np.zeros(psize)
    for j in range(psize):
      fitnesses[j] = pop[j].fitness_calc(problem)

  bestfitness = np.max(fitnesses)
  fitnesses = fitnesses
  fitness_hist[0, :] = fitnesses

  fitness_time[0] =  time() - start
  counter = 0 # Generation Repetition without fitness improvement counter
  for i in range(1, gen):
      start = time()

      # Crossover & Mutation
      offsprings = total_crossover(pop)
      for j in range(psize): offsprings[j].mutation(1/dims)

      # Fitness Calculation

      if multi_proc:
        with concurrent.futures.ProcessPoolExecutor(max_workers=workers) as executor:
          processes = []
          for j in range(psize):

              proc = executor.submit(offsprings[j].fitness_calc, problem)
              processes.append(proc)

          cfitnesses = [proc.result() for proc in processes]

          for j in range(psize):
              offsprings[j].fitness = cfitnesses[j]
      else:
        cfitnesses = np.zeros(psize)
        for j in range(psize):
          cfitnesses[j] = offsprings[j].fitness_calc(problem)

      # Selection
      pop, fitnesses = total_selection(np.concatenate((pop, offsprings)),
                                 np.concatenate((fitnesses, cfitnesses)), psize)

      fitness_hist[i, :] = fitnesses
      print('Generation %d best fitness = %f' % (i, bestfitness))

      if fitnesses[0] > bestfitness:
          bestfitness = fitnesses[0]
          counter = 0
      else:
          counter += 1
      fitness_time[i] = time() - start
      if counter == 20 and stop_condition:
          fitness_hist[i:, :] = fitnesses[0]
          break

  best_sol = pop[0]
  if create_model:
    model = ProbabilisticModel('umd')
    model.buildModel(Chromosome.genes_to_numpy(pop))
    src_models.append(model)

  return src_models, best_sol, fitness_hist, fitness_time

def continues_ea(fitness_func, init_func, dim, psize=100, gen=100, muc=10,
                           mum=10, stop_condition=True, create_model=True):

    src_model = None

    fitness_hist = np.zeros((gen, psize))
    fitness_time = np.zeros((gen))

    class ChromosomeEA(ChromosomePole):
        def fitness_calc(self):
            if self.fitness != float('-inf'):
                return self.fitness
            self.fitness = fitness_func(self.genes)
            return self.fitness

    pop = get_pop_init(psize, dim, init_func, p_type=ChromosomeEA)
    start = time()

    for j in range(psize):
        pop[j].fitness_calc()

    bestfitness = np.max(pop).fitness
    fitness = Chromosome.fitness_to_numpy(pop)
    fitness_hist[0, :] = fitness

    fitness_time[0] =  time() - start
    print('Generation %d best fitness = %f' % (0, bestfitness))
    counter = 0 # Generation Repetition without fitness improvement counter
    for i in range(1, gen):
        start = time()
        randlist = np.random.permutation(psize)
        offsprings = np.ndarray(psize, dtype=object)

        # Crossover & Mutation
        for j in range(0, psize, 2):
            offsprings[j] = ChromosomeEA(dim)
            offsprings[j+1] = ChromosomeEA(dim)
            p1 = randlist[j]
            p2 = randlist[j+1]
            offsprings[j].genes, offsprings[j+1].genes = sbx_crossover(pop[p1], pop[p2], muc, dim)
            offsprings[j].mutation(mum, dim)
            offsprings[j+1].mutation(mum, dim)

        # Fitness Calculation
        cfitness = np.zeros(psize)
        for j in range(psize):
            cfitness[j] = offsprings[j].fitness_calc()

        # Selection
        pop, fitness = total_selection(np.concatenate((pop, offsprings)),
                                 np.concatenate((fitness, cfitness)), psize)

        fitness_hist[i, :] = fitness

        if fitness[0] > bestfitness:
            bestfitness = fitness[0]
            counter = 0
        else:
            counter += 1

        print('Generation %d best fitness = %f' % (i, bestfitness))

        fitness_time[i] = time() - start

    best_sol = pop[0]
    if create_model:
        model = ProbabilisticModel('mvarnorm')
        print('build model input shape: ', ChromosomePole.genes_to_numpy(pop).shape)
        model.buildModel(ChromosomePole.genes_to_numpy(pop))
        print('model mean: ', model.mean)
        print("Model built successfully!")
        src_model = model
    elif not create_model:
        print("Evolutionary algorithm didn't reach the criteria!")

    return src_model, best_sol, fitness_hist, fitness_time

def transfer_continues_ea(fitness_func, init_func, dim, src_models, psize=100, gen=100,
                 muc=10, mum=10, reps=1, delta=2, build_model=True):

    if not src_models:
        raise ValueError('No probabilistic models stored for transfer optimization.')

    class ChromosomeEA(ChromosomePole):
        def fitness_calc(self):
            if self.fitness != float('-inf'):
                return self.fitness
            self.fitness = fitness_func(self.genes)
            return self.fitness

    fitness_hist = np.zeros([reps, gen, psize])
    fitness_time = np.zeros((reps, gen,))
    alpha = list()

    pop = None
    for rep in range(reps):
        alpha_rep = []
        pop = get_pop_init(psize, dim, init_func, p_type=ChromosomeEA)
        start = time()
        for j in range(psize):
            pop[j].fitness_calc()

        bestfitness = np.max(pop).fitness
        fitness = Chromosome.fitness_to_numpy(pop)
        fitness_hist[rep, 0, :] = fitness
        fitness_time[rep, 0] = time() - start
        print('Generation 0 best fitness = %f' % bestfitness)

        for i in range(1, gen):
            start = time()
            if i % delta == 0:
                mixModel = MixtureModel(src_models)
                mixModel.createTable(Chromosome.genes_to_numpy(pop), True, 'mvarnorm')
                mixModel.EMstacking()
                mixModel.mutate()
                offsprings = mixModel.sample(psize)
                offsprings = np.array([ChromosomeEA(offspring) for offspring in offsprings])
                alpha_rep = np.concatenate((alpha_rep, mixModel.alpha), axis=0)
#                 print('Mixture coefficients: %s' % np.array(mixModel.alpha))
            else:
                # Crossover & Mutation
                randlist = np.random.permutation(psize)
                offsprings = np.ndarray(psize, dtype=object)
                for j in range(0, psize, 2):
                    offsprings[j] = ChromosomeEA(dim)
                    offsprings[j+1] = ChromosomeEA(dim)
                    p1 = randlist[j]
                    p2 = randlist[j+1]
                    offsprings[j].genes, offsprings[j+1].genes = sbx_crossover(pop[p1], pop[p2], muc, dim)
                    offsprings[j].mutation(mum, dim)
                    offsprings[j+1].mutation(mum, dim)

            # Fitness Calculation
            cfitness = np.zeros(psize)
            for j in range(psize):
                cfitness[j] = offsprings[j].fitness_calc()

            if i % delta == 0:
                print('cfitness mean: ', np.mean(cfitness))
                print('cfitness max: ', np.max(cfitness))
                print('cfitness min: ', np.min(cfitness))

            # Selection
            pop, fitness = total_selection(np.concatenate((pop, offsprings)),
                                    np.concatenate((fitness, cfitness)), psize)

            fitness_hist[rep, i, :] = fitness
            fitness_time[rep, i] = time() - start

            if fitness[0] > bestfitness:
                bestfitness = fitness[0]

            print('Generation %d best fitness = %f' % (i, bestfitness))
            print('Generation %d mean fitness = %f' % (i, np.mean(fitness)))

        print()

        alpha.append(alpha_rep)

    model = None
    print('fitness_hist: ', fitness_hist[0, -1, 0])
    if build_model:
        model = ProbabilisticModel('mvarnorm')
        print('build model input shape: ', Chromosome.genes_to_numpy(pop).shape)
        model.buildModel(Chromosome.genes_to_numpy(pop))
        print('model mean: ', model.mean)
        print("Model built successfully!")
    else:
        print("Evolutionary algorithm didn't reach the criteria!")

    if build_model:
        return fitness_hist[0, ...], alpha, fitness_time[0, ...], model, np.max(pop).genes
    else:
        return fitness_hist, alpha, fitness_time

def transfer_ea(problem, dims, delta, psize=100, gen=100,
               create_model=True, stop_condition=True,
               src_models=[]):
  # load probabilistic models

  if src_models is None or len(src_models) == 0:
      raise ValueError('No probabilistic models stored for transfer optimization.')

  init_func = lambda n: np.round(np.random.rand(n))
  fitness_hist = np.zeros([gen, psize])
  fitness_time = np.zeros((gen))

  alpha_rep = []
  counter = 0

  pop = get_pop_init(psize, dims, init_func)
  start = time()
  for i in range(psize): pop[i].fitness_calc(problem)

  bestfitness = np.max(pop).fitness
  fitness = Chromosome.fitness_to_numpy(pop)
  fitness_hist[0, :] = fitness
  fitness_time[0] = time() - start
  print('Generation 0 best fitness = %f' % bestfitness)
  for i in range(1, gen):

      start = time()
      if i % delta == 0:
          mixModel = MixtureModel(src_models)
          mixModel.createTable(Chromosome.genes_to_numpy(pop), True, 'umd')
          mixModel.EMstacking()
          mixModel.mutate()
          offsprings = mixModel.sample(psize)

          offsprings = np.array([Chromosome(offspring) for offspring in offsprings])
          alpha_rep.append(mixModel.alpha)
          print('Mixture coefficients: %s' % np.array(mixModel.alpha))

      else:
          # Crossover & Mutation
          offsprings = total_crossover(pop)
          for j in range(psize): offsprings[j].mutation(1/dims)

      # Fitness Calculation
      cfitness = np.zeros(psize)
      for j in range(psize):
        cfitness[j] = offsprings[j].fitness_calc(problem)

      # Selection
      pop, fitness = total_selection(np.concatenate((pop, offsprings)),
                                np.concatenate((fitness, cfitness)), psize)

      fitness_hist[i, :] = fitness

      if fitness[0] > bestfitness:
          bestfitness = fitness[0]
          counter = 0
      else:
          counter += 1

      fitness_time[i] = time() - start

      if counter == 20 and stop_condition:
          fitness_hist[i:, :] = fitness[0]
          break

      print('Generation %d best fitness = %f' % (i, bestfitness))

  best_sol = pop[0]
  src_model = None
  if create_model:
    src_model = ProbabilisticModel('umd')
    print('build model input shape: ', Chromosome.genes_to_numpy(pop).shape)
    src_model.buildModel(Chromosome.genes_to_numpy(pop))
    print('probOne_noisy: ', src_model.probOne_noisy)
    print('probZero_noisy: ', src_model.probZero_noisy)

  return src_model, best_sol, fitness_hist, fitness_time

def evolutionary_algorithm_pole(sLen, psize=100, gen=100, muc=10,
                             mum=10, stop_condition=True, create_model=True,
                             multi_proc=False, workers=1):

  src_model = None

  fitness_hist = np.zeros((gen, psize))
  fitness_time = np.zeros((gen))

  cart = PoledCart(sLen)

  n_input = 6
  n_hidden = 10
  n_output = 1
  net = Net(n_input, n_hidden, n_output)
  n_vars = net.nVariables

  init_func = lambda n: 12 * np.random.rand(n) - 6
  pop = get_pop_init(psize, n_vars, init_func, p_type='double_pole')
  start = time()

  if multi_proc:
    with concurrent.futures.ProcessPoolExecutor(max_workers=workers) as executor:
      processes = []
      for j in range(psize):

          proc = executor.submit(pop[j].fitness_calc, net, cart, sLen)
          processes.append(proc)

      fitnesses = [proc.result() for proc in processes]

      for j in range(psize):
          pop[j].fitness = fitnesses[j]
  else:
    fitnesses = np.zeros(psize)
    for j in range(psize):
      fitnesses[j] = pop[j].fitness_calc(net, cart, sLen)

  bestfitness = np.max(fitnesses)
  fitness_hist[0, :] = fitnesses

  fitness_time[0] = start - time()
  counter = 0 # Generation Repetition without fitness improvement counter
  for i in range(1, gen):
      start = time()
      randlist = np.random.permutation(psize)
      offsprings = np.ndarray(psize, dtype=object)

      # Crossover & Mutation
      for j in range(0, psize, 2):
            offsprings[j] = ChromosomePole(n_vars)
            offsprings[j+1] = ChromosomePole(n_vars)
            p1 = randlist[j]
            p2 = randlist[j+1]
            offsprings[j].genes, offsprings[j+1].genes = sbx_crossover(pop[p1], pop[p2], muc, n_vars)
            offsprings[j].mutation(mum, n_vars)
            offsprings[j+1].mutation(mum, n_vars)

      # Fitness Calculation
      if multi_proc:
        with concurrent.futures.ProcessPoolExecutor(max_workers=workers) as executor:
          processes = []
          for j in range(psize):

              proc = executor.submit(offsprings[j].fitness_calc, net, cart, sLen)
              processes.append(proc)

          cfitnesses = [proc.result() for proc in processes]

          for j in range(psize):
              offsprings[j].fitness = cfitnesses[j]
      else:
        cfitnesses = np.zeros(psize)
        for j in range(psize):
          cfitnesses[j] = offsprings[j].fitness_calc(net, cart, sLen)

      cfitnesses = np.zeros(psize)
      for j in range(psize):
          # print(pop[j].genes)
          cfitnesses[j] = offsprings[j].fitness_calc(net, cart, sLen)

      # Selection
      pop, fitnesses = total_selection(np.concatenate((pop, offsprings)),
                                 np.concatenate((fitnesses, cfitnesses)), psize)

      fitness_hist[i, :] = fitnesses

      if fitnesses[0] > bestfitness:
          bestfitness = fitnesses[0]
          counter = 0
      else:
          counter += 1

      print('Generation %d best fitness = %f' % (i, bestfitness))
      if fitnesses[0] - 2000 > -0.0001 and stop_condition:
          print('Solution found!')
          fitness_hist[i:, :] = fitnesses[0]
          break

      fitness_time[i] = time() - start

  best_sol = pop[0]
  if create_model and fitness_hist[-1, 0] - 2000 > -0.0001:
    model = ProbabilisticModel('mvarnorm')
    print('build model input shape: ', Chromosome.genes_to_numpy(pop).shape)
    model.buildModel(Chromosome.genes_to_numpy(pop))
    print("Model built successfully!")
    src_model = model
  elif not create_model:
    print("Evolutionary algorithm didn't reach the criteria!")
    # src_models.append(model)

  return src_model, best_sol, fitness_hist, fitness_time

max_stars_repo_name: xlouba/Transfer-Optimization

max_issues_repo_name: xlouba/Transfer-Optimization

max_forks_repo_name: xlouba/Transfer-Optimization


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

import numpy as np
import cv2
img = cv2.imread('messi5.jpg',0)
rows,cols = img.shape

M = np.float32([[1,0,100],[0,1,50]])
dst = cv2.warpAffine(img,M,(cols,rows))

cv2.imshow('img',dst)
cv2.waitKey(0)
cv2.destroyAllWindows()

# ROTATION

img = cv2.imread('messi5.jpg',0)
rows,cols = img.shape

M = cv2.getRotationMatrix2D((cols/2,rows/2),90,1)
dst = cv2.warpAffine(img,M,(cols,rows))
cv2.imshow('img',dst)
cv2.waitKey(0)
cv2.destroyAllWindows()

max_stars_repo_path: OpenCV Python/4. Image Processing/3. geometric transformation/translation.py


text:
import skrf as rf
import numpy as np
import os
import xlsxwriter

def write_s1p(filename, distance, data):
    f = open('Files_Time/' + filename + ".s1p","w+")
    f.write('!Keysight Technologies N9923A: A.08.19\n')
    f.write('!Date: Thursday, 16 January 2020 11:30:00\n')
    f.write('!TimeZone: (GMT-03:00) Brasilia!Model: N9923A\n')
    f.write('!Serial: MY51491501\n')
    f.write('!GPS Latitude: \n')
    f.write('!GPS Longitude: \n')
    f.write('!GPS TimeStamp: 0001-01-01 00:00:00Z\n')
    f.write('!GPS Seconds Since Last Read: 0\n')
    f.write('!CHECKSUM:1078980270\n')
    f.write('!Correction: S11(ON U)\n')
    f.write('!S1P File: Measurement: S11:\n')
    f.write('# D S DB R 50\n')

    for i in range(len(distance)):
        f.write(str(distance[i]) + ' ' + str(data[i]) + ' ' + str(0)  + '\n')

    f.close()

def smooth(y, box_pts):
    box = np.ones(box_pts)/box_pts
    y_smooth = np.convolve(y, box, mode='same')
    return y_smooth

wb = xlsxwriter.Workbook('Time01.xlsx')
ws = wb.add_worksheet()
head_writed = False

for i,file1 in enumerate(os.listdir("Files")):

    if file1.endswith(".s1p"):
        s = file1.split('_')
        index1 = s[1]
        fonte = s[2]
        haste = s[3][0:-4]
        f = open('Files_Time/' + file1 + ".s1p","r+")
        s = f.readline()
        j = 0
        while s:
            if s.startswith('!') or s.startswith('#'):
                pass
            else:
                l = s.split(' ')
                if head_writed is False:
                    ws.write(i,j,np.round(float(l[0]),2))
                    ws.write(i+1,j,l[1])
                else:
                    ws.write(i+1,j,l[1])
                j = j + 1
            s = f.readline()
        head_writed = True
    #df.to_excel('Modulo de S11.xlsx', sheet_name='sheet1', index=False)
wb.close()
##for i,file1 in enumerate(os.listdir("Files")):
##
##    if file1.endswith(".s1p"):
##        probe = rf.Network('Files/'+ file1)
##        x,y = probe.plot_s_db_time(window=('kaiser', 6))
##        x = 0.92*3*10**8*x/2
##        y = smooth(y,3)
##        write_s1p(file1, x, y )


text:
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
import pysindy as ps

# ignore user warnings
import warnings
warnings.filterwarnings("ignore", category=UserWarning)

np.random.seed(42)  # Seed for reproducibility

# Linear, damped harmonic oscillator:
def linear_damped_SHO(t, x):
	"""x' = -0.1x+2y
	y' = -2x-0.1y
	"""
	return [-0.1 * x[0] + 2 * x[1], -2 * x[0] - 0.1 * x[1]]

# Generate training data
delta_t = 0.01
t_train = np.arange(0, 20, delta_t)
t_train_range = (t_train[0], t_train[-1])
x0_train = [1, 0] # initial condition
x_train = solve_ivp(linear_damped_SHO, t_train_range, x0_train, t_eval=t_train, rtol = 1e-12, method = 'LSODA', atol = 1e-12).y.T

poly_order = 5 # looking at space of polynomials up to 5th order
threshold = 0.05 # selecting by default a threshold equal t0 0.5

# Build the model
model = ps.SINDy(
    optimizer=ps.STLSQ(threshold=threshold),
    feature_library=ps.PolynomialLibrary(degree=poly_order),
)

# Look at the results over multiple noise levels
x_dot_train_measured = np.array(
    [linear_damped_SHO(0, x_train[i]) for i in range(t_train.size)]
)
noise_levels = [1e-4, 1e-3, 1e-2, 1e-1, 1]
models = []
t_sim = np.arange(0, 20, delta_t)
x_sim = []

for eps in noise_levels:
    model = ps.SINDy(
        optimizer=ps.STLSQ(threshold=threshold),
        feature_library=ps.PolynomialLibrary(degree=poly_order),
    )
    model.fit(
        x_train,
        t=delta_t,
        x_dot=x_dot_train_measured
        + np.random.normal(scale=eps, size=x_train.shape),
        quiet=True,
    )
    models.append(model)
    x_sim.append(model.simulate(x_train[0], t_sim))

# show model with lowest noise level:
print(f'Identified model with noise level {noise_levels[0]} is:')
models[0].print()

# Look at the error over multiple thresholds
thresholds = [i/100 for i in range(1,100,5)]
models_t = []
t_sim = np.arange(0, 20, delta_t)
x_simt = []
eps = 0.001
for t in thresholds:
    model = ps.SINDy(
        optimizer=ps.STLSQ(threshold=t),
        feature_library=ps.PolynomialLibrary(degree=poly_order),
    )
    model.fit(
        x_train,
        t=delta_t,
        x_dot=x_dot_train_measured
        + np.random.normal(scale=eps, size=x_train.shape),
        quiet=True,
    )
    models_t.append(model)
    x_simt.append(model.simulate(x_train[0], t_sim))

# Plot results

# plot actual training data without noise
fig = plt.figure(figsize=(15, 4))
ax = fig.add_subplot(131)
ax.plot(
    x_train[: t_sim.size, 0],
    x_train[: t_sim.size, 1]
)
plt.title("Training data no noise")
ax.set(xlabel="$x$", ylabel="$y$")

# plot identified model data with small noise
model_idx = 0
ax = fig.add_subplot(132)
ax.plot(x_sim[model_idx][:, 0], x_sim[model_idx][:, 1])
plt.title(f"SINDy identified system, Noise level={noise_levels[model_idx]}")
ax.set(xlabel="$x$", ylabel="$y$")

# plot identified data with higher noise
model_idx = 3
ax = fig.add_subplot(133)
ax.plot(x_sim[model_idx][:, 0], x_sim[model_idx][:, 1])
plt.title(f"SINDy identified system, Noise level={noise_levels[model_idx]}")
ax.set(xlabel="$x$", ylabel="$y$")
plt.savefig('SINDy/graphs/vibration.png')

fig = plt.figure(figsize=(12, 5))
model_idx = 0
ax = fig.add_subplot(221)
ax.plot(t_sim, x_train[: t_sim.size, 0], "r", label = 'Training data')
ax.plot(t_sim, x_sim[model_idx][:, 0], "b--", label = 'Simulated data')
plt.title(f"x, y vs t for Noise level={noise_levels[model_idx]}")
plt.ylabel("x")
plt.legend()

ax = fig.add_subplot(223)
ax.plot(t_sim, x_train[: t_sim.size, 1], "r",label = 'Training data')
ax.plot(t_sim, x_sim[model_idx][:, 1], "b--", label = 'Simulated data')
plt.xlabel("time")
plt.ylabel("y")
plt.legend()

model_idx = 3
ax = fig.add_subplot(222)
ax.plot(t_sim, x_train[: t_sim.size, 0], "r",label = 'Training data')
ax.plot(t_sim, x_sim[model_idx][:, 0], "b--", label = 'Simulated data')
plt.title(f"x, y vs t for Noise level={noise_levels[model_idx]}")
plt.ylabel("x")
plt.legend()

ax = fig.add_subplot(224)
ax.plot(t_sim, x_train[: t_sim.size, 1], "r",label = 'Training data')
ax.plot(t_sim, x_sim[model_idx][:, 1], "b--", label = 'Simulated data')
plt.xlabel("time")
plt.ylabel("y")
plt.legend()
plt.savefig('SINDy/graphs/vibration2.png')

# errors (RMSE vs noise):
rmse_average_noise = []
for model_idx in range(len(noise_levels)):
	rmse_average_noise.append(np.average(np.sqrt((x_sim[model_idx][:, 0]-x_train[: t_sim.size, 0])**2+(x_sim[model_idx][:, 1]-x_train[: t_sim.size, 1])**2),axis = 0))
fig,ax = plt.subplots(figsize=(12, 8))
ax.plot(noise_levels, rmse_average_noise, '-ok')
ax.set_title('RMSE vs noise')
ax.set_ylabel('RMSE in both x and y')
ax.set_xlabel('Noise levels')
ax.set_yscale('log')
plt.savefig('SINDy/graphs/vibration3.png')
# fig.tight_layout()

# errors (RMSE vs threshold):
rmse_average_t = []
for model_idx in range(len(thresholds)):
	rmse_average_t.append(np.average(np.sqrt((x_simt[model_idx][:, 0]-x_train[: t_sim.size, 0])**2+(x_simt[model_idx][:, 1]-x_train[: t_sim.size, 1])**2),axis = 0))
fig,ax = plt.subplots(figsize=(12, 8))
ax.plot(thresholds, rmse_average_t, '-ok')
ax.set_title('RMSE vs threshold level')
ax.set_ylabel('RMSE in both x and y')
ax.set_xlabel('Threshold levels')
ax.set_yscale('log')
plt.savefig('SINDy/graphs/vibration4.png')

# plt.show()

reason: import numpy,from scipy


text:
import asyncio
import math
import time
import wave

import cv2
import numpy

from ..mediastreams import (AudioFrame, AudioStreamTrack, VideoFrame,
                            VideoStreamTrack)

AUDIO_PTIME = 0.020  # 20ms audio packetization

def frame_from_bgr(data_bgr):
    data_yuv = cv2.cvtColor(data_bgr, cv2.COLOR_BGR2YUV_I420)
    return VideoFrame(width=data_bgr.shape[1], height=data_bgr.shape[0], data=data_yuv.tobytes())

def frame_from_gray(data_gray):
    data_bgr = cv2.cvtColor(data_gray, cv2.COLOR_GRAY2BGR)
    data_yuv = cv2.cvtColor(data_bgr, cv2.COLOR_BGR2YUV_I420)
    return VideoFrame(width=data_bgr.shape[1], height=data_bgr.shape[0], data=data_yuv.tobytes())

def frame_to_bgr(frame):
    data_flat = numpy.frombuffer(frame.data, numpy.uint8)
    data_yuv = data_flat.reshape((math.ceil(frame.height * 12 / 8), frame.width))
    return cv2.cvtColor(data_yuv, cv2.COLOR_YUV2BGR_I420)

class AudioFileTrack(AudioStreamTrack):
    """
    An AudioStreamTrack subclass for reading audio from a WAV file.
    """
    def __init__(self, path):
        self.last = None
        self.reader = wave.open(path, 'rb')
        self.frames_per_packet = int(self.reader.getframerate() * AUDIO_PTIME)

    async def recv(self):
        # as we are reading audio from a file and not using a "live" source,
        # we need to control the rate at which audio is sent
        if self.last:
            now = time.time()
            await asyncio.sleep(self.last + AUDIO_PTIME - now)
        self.last = time.time()

        return AudioFrame(
            channels=self.reader.getnchannels(),
            data=self.reader.readframes(self.frames_per_packet),
            sample_rate=self.reader.getframerate())

class VideoFileTrack(VideoStreamTrack):
    """
    A VideoStreamTrack subclass for reading video from a file.
    """
    def __init__(self, path):
        self.cap = cv2.VideoCapture(path)
        self.last = None
        self.ptime = 1 / self.cap.get(cv2.CAP_PROP_FPS)

    async def recv(self):
        # as we are reading audio from a file and not using a "live" source,
        # we need to control the rate at which audio is sent
        if self.last:
            now = time.time()
            await asyncio.sleep(self.last + self.ptime - now)
        self.last = time.time()

        ret, frame = self.cap.read()
        if not ret:
            # loop
            self.cap.set(cv2.CAP_PROP_POS_FRAMES, 0)
            ret, frame = self.cap.read()

        return frame_from_bgr(frame)

max_stars_repo_name: nickaknudson/aiortc

max_issues_repo_name: nickaknudson/aiortc

max_forks_repo_name: nickaknudson/aiortc


text:
import os

import numpy as np

exps = dict()
exps["1.a.1"] = dict(
    name="fd",
    commit="231423",
    exec="""
        fsdf
    """,
    haha="fad",
)

max_stars_repo_name: klauscc/vedatad

max_issues_repo_name: klauscc/vedatad

max_forks_repo_name: klauscc/vedatad


text:
"""
Code to work with ocnN.  It gets the data_dict entries for one time and
Saves the result.  Meant to work as a parallel subprocess
"""

from lo_tools import zfun
import xarray as xr
import numpy as np
import sys
from scipy.spatial import cKDTree
import argparse
import pickle

parser = argparse.ArgumentParser()
parser.add_argument('-grid_fn', type=str) # full path to nest grid.nc file
parser.add_argument('-his_fn', type=str) # full path to history file to interpolate from
parser.add_argument('-start_type', type=str) # continuation or new
parser.add_argument('-out_fn', type=str) # full path to output pickle file
args = parser.parse_args()

def get_bounds(x_big, y_big, x_small, y_small, pad=3):
    """
    This function takes two pairs of plaid, 2-D, lon, lat arrays:
    - one pair bigger (that we hope to nest inside) and
    - one pair smaller (the grid of the nest)
    and returns the indices to use for making a trimmed version of
    the bigger grid that the smaller grid still fits inside.
    We add "pad" around the edges to make sure things fit comfortably.
    """
    # First: error checking
    if (x_small[0,0] < x_big[0,pad]) or (x_small[0,-1] > x_big[0,-pad]):
        print('ERROR: lon out of bounds ')
        sys.exit()
    if (y_small[0,0] < y_big[pad,0]) or (y_small[-1,0] > y_big[-pad,0]):
        print('ERROR: lat out of bounds ')
        sys.exit()
    # Second: get indices
    ix0 = zfun.find_nearest_ind(x_big[0,:], x_small[0,0]) - pad
    ix1 = zfun.find_nearest_ind(x_big[0,:], x_small[0,-1]) + pad
    iy0 = zfun.find_nearest_ind(y_big[:,0], y_small[0,0]) - pad
    iy1 = zfun.find_nearest_ind(y_big[:,0], y_small[-1,0]) + pad
    return ix0, ix1, iy0, iy1

tag_list = ['rho', 'u', 'v']

# the new grid
ds = xr.open_dataset(args.grid_fn)

xx = {}; yy = {}; mm = {}; xynew = {}
for tag in tag_list:
    xx[tag] = ds['lon_' + tag].values
    yy[tag] = ds['lat_' + tag].values
    mm[tag] = ds['mask_' + tag].values

ds.close()

if args.start_type == 'continuation':
    pad = 20
elif args.start_type == 'new':
    pad = 0
else:
    print('Error: Unrecognized start_type')
    sys.exit()
if pad > 0:
    # mask out the inside of the nest fields, since we only use
    # the edges (unless Ldir['start_type']=='new')
    for tag in tag_list:
        mm[tag][pad:-pad, pad:-pad] = 0 # this speeds things up

for tag in tag_list:
    xynew[tag] = np.array((xx[tag][mm[tag]==1],yy[tag][mm[tag]==1])).T

# Create 2-D search trees for the old grid
ds = xr.open_dataset(args.his_fn)
N = len(ds.s_rho.values)
xtrim = {}; ytrim = {}; mtrim = {}; xyT = {}
ix0 = {}; ix1 = {}; iy0 = {}; iy1 = {}
for tag in tag_list:
    x = ds['lon_' + tag].values
    y = ds['lat_' + tag].values
    m = ds['mask_' + tag].values # 1=water
    # trim the old grid before making the search tree
    ix0[tag], ix1[tag], iy0[tag], iy1[tag] = get_bounds(x, y, xx[tag], yy[tag])
    xtrim[tag] = x[iy0[tag]:iy1[tag], ix0[tag]:ix1[tag]]
    ytrim[tag] = y[iy0[tag]:iy1[tag], ix0[tag]:ix1[tag]]
    mtrim[tag] = m[iy0[tag]:iy1[tag], ix0[tag]:ix1[tag]]
    xyorig = np.array((xtrim[tag][mtrim[tag]==1],ytrim[tag][mtrim[tag]==1])).T
    xyT[tag] = cKDTree(xyorig)
ds.close()

# associate variables to process with grids
vn_dict = {'salt':('rho',3), 'temp':('rho',3), 'zeta':('rho',2),
        'u':('u',3), 'v':('v',3), 'ubar':('u',2), 'vbar':('v',2)}

# create blank arrays for results
data_dict = dict()
for vn in vn_dict.keys():
    tag = vn_dict[vn][0]
    dm = vn_dict[vn][1]
    NR, NC = xx[tag].shape
    if dm == 2:
        data_dict[vn] = np.nan* np.ones((NR,NC))
    elif dm == 3:
        data_dict[vn] = np.nan* np.ones((N,NR,NC))

# Interpolate to fill all data arrays for new grid.
ds = xr.open_dataset(args.his_fn, decode_times=False)
data_dict['ocean_time'] = ds.ocean_time.values[0]
for vn in vn_dict.keys():
    tag = vn_dict[vn][0]
    dm = vn_dict[vn][1]
    if dm == 2:
        vtrim = ds[vn][0,iy0[tag]:iy1[tag], ix0[tag]:ix1[tag]].values
        vv = np.nan * np.ones(xx[tag].shape)
        vv[mm[tag]==1] = vtrim[mtrim[tag]==1][xyT[tag].query(xynew[tag], workers=-1)[1]]
        # note that "workers" has replaced "n_jobs"
        data_dict[vn][:, :] = vv
    elif dm == 3:
        for nn in range(N):
            vtrim = ds[vn][0,nn,iy0[tag]:iy1[tag], ix0[tag]:ix1[tag]].values
            vv = np.nan * np.ones(xx[tag].shape)
            vv[mm[tag]==1] = vtrim[mtrim[tag]==1][xyT[tag].query(xynew[tag], workers=-1)[1]]
            data_dict[vn][nn, :, :] = vv
ds.close()

pickle.dump(data_dict, open(args.out_fn,'wb'))

max_stars_repo_name: parkermac/LO

max_issues_repo_name: parkermac/LO

max_forks_repo_name: parkermac/LO


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

"""

Convnet training/fine-tuning.

Created on Thu Jul  4 23:37:36 2019

@author: vlado
"""
import os
import pickle
import random
import sklearn
import argparse
import numpy as np
import tensorflow as tf
from tensorflow.keras import Model
import tensorflow.keras.backend as K
from tensorflow.keras import optimizers
from tensorflow.keras.applications.resnet50 import ResNet50
from tensorflow.keras.callbacks import Callback, ModelCheckpoint
from tensorflow.keras.layers import Dense, GlobalAveragePooling2D

os.environ['PYTHONHASHSEED'] = '42'
np.random.seed(42)
random.seed(42)
tf.random.set_seed(42)

class LRScheduleSteps(Callback):
  def __init__(self, opt, lr, w, p, factor=10.):
    super(LRScheduleSteps, self).__init__()
    self.optimizer = opt
    self.lr = lr
    self.w = w
    self.p = p
    self.factor = factor
    self.history = {}

    self.step = 0

  def on_train_begin(self, logs={}):
      logs = logs or {}

  def on_batch_end(self, batch, logs=None):
    logs = logs or {}

    if self.step < self.w:
      lr = self.lr / self.w * self.step
    else:
      lr = K.get_value(self.optimizer.lr)
      for step in self.p:
        if step == self.step:
          lr = lr / self.factor

    self.history.setdefault('lr', []).append(K.get_value(self.optimizer.lr))
    self.history.setdefault('iterations', []).append(self.step)

    for k, v in logs.items():
        self.history.setdefault(k, []).append(v)

    K.set_value(self.optimizer.lr, lr)
    self.step += 1

def preprocess(image, preproc):
    if preproc:
      image = tf.cast(image, tf.float32)
      means = tf.constant(np.reshape([123.68, 116.779, 103.939], (1, 1, 3)),
                          dtype=tf.float32)
      image = tf.math.subtract(image, means)
    else:
      image = tf.image.convert_image_dtype(image, tf.float32)

    return image

class LoadPreprocessImage():
  def __init__(self, image_path,
               load_size=(256, 256),
               dim=(224, 224),
               crop_size=(32, 256),
               preproc=False):
    self.image_path = image_path
    self.load_size = load_size
    self.dim = dim
    self.crop_size = crop_size
    self.preproc = preproc

  def __call__(self, record):
    image = tf.io.read_file(record['filename'])
    image = tf.image.decode_jpeg(image)
    image = preprocess(image, self.preproc)
    image = tf.image.resize(image, self.load_size)
    image = tf.image.random_crop(image, self.dim+(3,))
    n = tf.random.uniform((), maxval=4, dtype=tf.int32)
    image = tf.image.rot90(image, k=n)
    if tf.random.uniform(()) > 0.5:
      image = tf.image.flip_left_right(image)
    if tf.random.uniform(()) > 0.5:
      image = tf.image.flip_up_down(image)

    return image, record['class']

class LoadPreprocessImageVal():
  def __init__(self, image_path,
               load_size=(256, 256),
               dim=(224, 224),
               preproc=False):
    self.image_path = image_path
    self.load_size = load_size
    self.dim = dim
    self.preproc = preproc

  def __call__(self, record):
    image = tf.io.read_file(record['filename'])
    image = tf.image.decode_jpeg(image)
    image = preprocess(image, self.preproc)
    image = tf.image.convert_image_dtype(image, tf.float32)
    image = tf.image.resize(image, self.load_size)
    image = tf.image.central_crop(image, self.dim[0]/self.load_size[0])

    return image, record['class']

parser = argparse.ArgumentParser()
parser.add_argument('-p', '--path', dest='path', required=True, help='Path to the images')
parser.add_argument('-d', '--data-split', dest='data', required=True, help='Dataset split')
parser.add_argument('-n', '--name', dest='name', help='Model name')
parser.add_argument('-m', '--model', dest='model', required=False,
                    help='None/imagenet/path_to_the_pre-trained_model.')
parser.add_argument('--lr', dest='max_lr', type=float, required=True, help='Maximal learning rate.')
parser.add_argument('-b', '--batch-size', dest='batch_size', type=int,
                    required=False, default=100,
                    help='Batch size.')
args = parser.parse_args()

dataset = args.data
batch_size = args.batch_size
preproc = args.model is not None and 'imagenet' in args.model

with open(f'data_splits/{dataset}-split.pkl', 'rb') as f:
    data_partition = pickle.load(f)
with open(f'data_splits/{dataset}-le.pkl', 'rb') as f:
    le = pickle.load(f)

data_partition['train']['filename'] = [os.path.join(args.path, fname) for fname in data_partition['train']['filename']]
data_partition['test']['filename'] = [os.path.join(args.path, fname) for fname in data_partition['test']['filename']]

strategy = tf.distribute.MirroredStrategy()

nr_classes = len(le.classes_)
nr_training = len(data_partition['train']['class'])
steps_per_epoch = nr_training // batch_size

options = tf.data.Options()
options.experimental_distribute.auto_shard_policy = tf.data.experimental.AutoShardPolicy.OFF

train_list_ds = tf.data.Dataset.from_tensor_slices(data_partition['train'])
train_ds = train_list_ds.shuffle(nr_training)
train_ds = train_ds.map(LoadPreprocessImage(args.path,
                                            load_size=(292, 292),
                                            dim=(256, 256),
                                            crop_size=(32, 256),
                                            preproc=preproc),
                        num_parallel_calls=tf.data.experimental.AUTOTUNE)
train_ds = train_ds.batch(batch_size)
train_ds = train_ds.prefetch(tf.data.experimental.AUTOTUNE)
train_ds = train_ds.with_options(options)

test_list_ds = tf.data.Dataset.from_tensor_slices(data_partition['test'])
test_ds = test_list_ds.map(LoadPreprocessImageVal(args.path,
                                                  load_size=(292, 292),
                                                  dim=(256, 256),
                                                  preproc=preproc),
                           num_parallel_calls=tf.data.experimental.AUTOTUNE)
test_ds = test_ds.batch(batch_size)
test_ds = test_ds.prefetch(tf.data.experimental.AUTOTUNE)
test_ds = test_ds.with_options(options)

with strategy.scope():
    if args.model is None or args.model == 'imagenet':
        base_model = ResNet50(include_top=False,
                              weights=args.model,
                              input_shape=(256, 256, 3))
        x = base_model.output
        x = GlobalAveragePooling2D()(x)
    else:
        base_model = tf.keras.models.load_model(args.model)
        x = base_model.layers[-1].input

    x = Dense(nr_classes, activation='softmax')(x)
    clf = Model(base_model.input, x)

    opt = optimizers.Adam(learning_rate=1e-5)
    clf.compile(opt,
                loss='sparse_categorical_crossentropy',
                metrics=['accuracy'])

lrsched = LRScheduleSteps(opt, args.max_lr, 5*steps_per_epoch,
                          [50*steps_per_epoch, 70*steps_per_epoch, 90*steps_per_epoch], factor=5.)

checkpoint = ModelCheckpoint('models/checkpoint.h5',
                             save_freq=10*batch_size,
                             save_best_only=False)

h = clf.fit(train_ds,
            epochs=100,
            callbacks=[checkpoint, lrsched],
            validation_data=test_ds)

filename = 'models/rssc_resnet50'
if args.model == 'imagenet':
    filename += '_imagenet'
if args.name is not None:
    filename += '_{}'.format(args.name)
filename += '_{}'.format(dataset)
if args.model is not None:
    filename += '_ft'
filename += '.h5'

clf.save(filename)

loss, acc = clf.evaluate(test_ds)

print('Classification accuracy: {:.2f}'.format(100*acc))

max_stars_repo_name: risojevicv/RSSC-transfer

max_issues_repo_name: risojevicv/RSSC-transfer

max_forks_repo_name: risojevicv/RSSC-transfer


text:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

def import_txt(filename):
    data = pd.read_csv(filename, sep="\t", header=None)

    data = data.sample(frac=1, random_state=42)

    X = data.iloc[:, 2:]
    y = data.iloc[:, 1]

    unique_labels = list(set(y))

    return np.array(X), np.array(y), unique_labels

# TODO: implement from scratch

def _make_incidence_matrix(labels):
    incidence_matrix = np.zeros((len(labels), len(labels)))

    for i in range(len(labels)):
        for j in range(len(labels)):
            incidence_matrix[i][j] = 1 if labels[i] == labels[j] else 0
    return incidence_matrix

def _get_confusion_matrix(gt_matrix, p_matrix):
    TT = TF = FT = FF = 0

    for i in range(len(p_matrix)):
        for j in range(len(gt_matrix)):

            if p_matrix[i][j] != gt_matrix[i][j]:
                if p_matrix[i][j] == 1:
                    TF += 1
                else:
                    FT += 1
            else:
                if p_matrix[i][j] == 1:
                    TT += 1
                else:
                    FF += 1

    return TT, TF, FT, FF

def jaccard(ground_truth, predicted):
    gt_matrix = _make_incidence_matrix(ground_truth)
    p_matrix = _make_incidence_matrix(predicted)

    TT, TF, FT, _ = _get_confusion_matrix(gt_matrix, p_matrix)
    return float(TT) / (TT + TF + FT)

def rand(ground_truth, predicted):
    gt_matrix = _make_incidence_matrix(ground_truth)
    p_matrix = _make_incidence_matrix(predicted)

    TT, TF, FT, FF = _get_confusion_matrix(gt_matrix, p_matrix)
    return float(TT + FF) / (TT + TF + FT + FF)

def pca(X):
    pca = PCA(n_components=2)
    pca.fit(X)
    new_X = pca.transform(X)

    return new_X

def scatter(X1, X2, labels, unique_labels, algo="", file_name=""):
    plt.figure(figsize=[17, 12])
    if algo and file_name:
        plt.title("Algorithm: {} Dataset: {}".format(algo, file_name))
    unique_encoded = [i for i in range(len(unique_labels))]

    colors = [plt.cm.jet(float(i)/max(unique_encoded))
              for i, u in enumerate(unique_labels)]

    for i, u in enumerate(unique_labels):
        xi = [X1[j] for j in range(len(X1)) if labels[j] == u]
        yi = [X2[j] for j in range(len(X2)) if labels[j] == u]
        plt.scatter(xi, yi, c=colors[i], label=str(u), s=35)

    plt.legend()

    plt.show()

max_stars_repo_name: socket-var/clustering-algorithms

max_issues_repo_name: socket-var/clustering-algorithms

max_forks_repo_name: socket-var/clustering-algorithms


text:
from sklearn.base import BaseEstimator
from dA import *
import numpy

import theano
import theano.tensor as T
from theano.tensor.shared_randomstreams import RandomStreams

class DenoisingAutoencoder(BaseEstimator):
    def __init__(self, n_hidden, learning_rate = 0.1, training_epochs = 20, corruption_level = 0.3, batch_size = 20, verbose=False):
        self.n_visible = None
        self.n_hidden = n_hidden
        self.learning_rate = learning_rate
        self.corruption_level = corruption_level
        self.training_epochs = training_epochs
        self.batch_size = batch_size
        self.verbose = verbose

        self.da = None
        self.x = T.matrix('x')

    def load_data(self,X):
        try:
            matrix = X.as_matrix()
        except AttributeError:
            matrix = X
        shared_x = theano.shared(numpy.asarray(matrix, dtype=theano.config.floatX), borrow=True)
        return shared_x

    def fit(self,X):
        self.n_visible = X.shape[1]

        train_set_x = self.load_data(X)

        # compute number of minibatches for training, validation and testing
        n_train_batches = train_set_x.get_value(borrow=True).shape[0] / self.batch_size

        # allocate symbolic variables for the data
        index = T.lscalar()    # index to a [mini]batch

        rng = numpy.random.RandomState(123)
        theano_rng = RandomStreams(rng.randint(2 ** 30))

        self.da = dA(
            numpy_rng=rng,
            theano_rng=theano_rng,
            input=self.x,
            n_visible=self.n_visible,
            n_hidden=self.n_hidden
        )

        cost, updates = self.da.get_cost_updates(
            corruption_level=self.corruption_level,
            learning_rate=self.learning_rate
        )

        train_da = theano.function(
            [index],
            cost,
            updates=updates,
            givens={
                self.x: train_set_x[index * self.batch_size: (index + 1) * self.batch_size]
            }
        )

        # training...
        for epoch in xrange(self.training_epochs):
            c = []
            for batch_index in xrange(n_train_batches):
                c.append(train_da(batch_index))

            if self.verbose:
                print 'Training epoch %d, cost ' % epoch, numpy.mean(c)

    def transform(self,X):
        z = self.da.get_prediction()
        predict_da = theano.function([self.x],z)
        return predict_da(X)

    def transform_latent_representation(self,X):
        h = self.da.get_latent_representation()
        predict_da = theano.function([self.x],h)
        return predict_da(X)

reason: import numpy,import theano,from theano


text:
import h5py
import numpy as np

from ef.field.solvers.field_solver import FieldSolver
from ef.util.serializable_h5 import SerializableH5

class Simulation(SerializableH5):

    def __init__(self, time_grid, spat_mesh, inner_regions,
                 particle_sources,
                 electric_fields, magnetic_fields, particle_interaction_model,
                 output_filename_prefix, outut_filename_suffix, max_id=-1, particle_arrays=()):
        self.time_grid = time_grid
        self.spat_mesh = spat_mesh
        self.inner_regions = inner_regions
        self._field_solver = FieldSolver(spat_mesh, inner_regions)
        self.particle_sources = particle_sources
        self.electric_fields = electric_fields
        self.magnetic_fields = magnetic_fields
        self.particle_interaction_model = particle_interaction_model
        self._output_filename_prefix = output_filename_prefix
        self._output_filename_suffix = outut_filename_suffix
        self.max_id = max_id
        self.particle_arrays = list(particle_arrays)

    @classmethod
    def init_from_h5(cls, h5file, filename_prefix, filename_suffix):
        domain = cls.load_h5(h5file)
        domain._output_filename_prefix = filename_prefix
        domain._output_filename_suffix = filename_suffix
        return domain

    def start_pic_simulation(self):
        self.eval_and_write_fields_without_particles()
        for src in self.particle_sources:
            particles = src.generate_initial_particles()
            if len(particles.ids):
                particles.ids = self.generate_particle_ids(len(particles.ids))
                self.particle_arrays.append(particles)
        self.prepare_recently_generated_particles_for_boris_integration()
        self.write_step_to_save()
        self.run_pic()

    def continue_pic_simulation(self):
        self.run_pic()

    def run_pic(self):
        total_time_iterations = self.time_grid.total_nodes - 1
        current_node = self.time_grid.current_node
        for i in range(current_node, total_time_iterations):
            print("Time step from {:d} to {:d} of {:d}".format(
                i, i + 1, total_time_iterations))
            self.advance_one_time_step()
            self.write_step_to_save()

    def prepare_recently_generated_particles_for_boris_integration(self):
        if self.particle_interaction_model.pic:
            self.eval_charge_density()
            self.eval_potential_and_fields()
        self.shift_new_particles_velocities_half_time_step_back()

    def advance_one_time_step(self):
        self.push_particles()
        self.apply_domain_constrains()
        if self.particle_interaction_model.pic:
            self.eval_charge_density()
            self.eval_potential_and_fields()
        self.update_time_grid()

    def eval_charge_density(self):
        self.spat_mesh.clear_old_density_values()
        self.spat_mesh.weight_particles_charge_to_mesh(self.particle_arrays)

    def eval_potential_and_fields(self):
        self._field_solver.eval_potential(self.spat_mesh, self.inner_regions)
        self._field_solver.eval_fields_from_potential(self.spat_mesh)

    def push_particles(self):
        self.boris_integration(self.time_grid.time_step_size)

    def apply_domain_constrains(self):
        # First generate then remove.
        # This allows for overlap of source and inner region.
        self.generate_new_particles()
        self.apply_domain_boundary_conditions()
        self.remove_particles_inside_inner_regions()

    def boris_integration(self, dt):
        for particles in self.particle_arrays:
            total_el_field, total_mgn_field = \
                self.compute_total_fields_at_positions(particles.positions)
            if total_mgn_field is not None and total_mgn_field.any():
                particles.boris_update_momentums(dt, total_el_field, total_mgn_field)
            else:
                particles.boris_update_momentum_no_mgn(dt, total_el_field)
            particles.update_positions(dt)

    def prepare_boris_integration(self, minus_half_dt):
        # todo: place newly generated particle_arrays into separate buffer
        for particles in self.particle_arrays:
            if not particles.momentum_is_half_time_step_shifted:
                total_el_field, total_mgn_field = \
                    self.compute_total_fields_at_positions(particles.positions)
                if total_mgn_field is not None and total_mgn_field.any():
                    particles.boris_update_momentums(minus_half_dt, total_el_field, total_mgn_field)
                else:
                    particles.boris_update_momentum_no_mgn(minus_half_dt, total_el_field)
                particles.momentum_is_half_time_step_shifted = True

    def compute_total_fields_at_positions(self, positions):
        total_el_field = np.zeros_like(positions)  # make sure shape is set, as += operators can't broadcast left side
        total_el_field += sum(f.get_at_points(positions, self.time_grid.current_time) for f in self.electric_fields)
        if self.particle_interaction_model.noninteracting:
            if self.inner_regions or not self.spat_mesh.is_potential_equal_on_boundaries():
                total_el_field += self.spat_mesh.field_at_position(positions)
        elif self.particle_interaction_model.binary:
            total_el_field += self.binary_electric_field_at_positions(positions)
            if self.inner_regions or not self.spat_mesh.is_potential_equal_on_boundaries():
                total_el_field += self.spat_mesh.field_at_position(positions)
        elif self.particle_interaction_model.pic:
            total_el_field += self.spat_mesh.field_at_position(positions)
        mgn_field = None
        if self.magnetic_fields:
            mgn_field = sum(f.get_at_points(positions, self.time_grid.current_time) for f in self.magnetic_fields)
        return total_el_field, mgn_field

    def binary_electric_field_at_positions(self, positions):
        return sum(
            np.nan_to_num(p.field_at_points(positions)) for p in self.particle_arrays)

    #
    # Push particles
    #

    def shift_new_particles_velocities_half_time_step_back(self):
        minus_half_dt = -1.0 * self.time_grid.time_step_size / 2.0
        #
        self.prepare_boris_integration(minus_half_dt)

    #
    # Apply domain constrains
    #

    def apply_domain_boundary_conditions(self):
        for arr in self.particle_arrays:
            collisions = self.out_of_bound(arr)
            arr.remove(collisions)
        self.particle_arrays = [a for a in self.particle_arrays if len(a.ids) > 0]

    def remove_particles_inside_inner_regions(self):
        for region in self.inner_regions:
            for p in self.particle_arrays:
                region.collide_with_particles(p)
            self.particle_arrays = [a for a in self.particle_arrays if len(a.ids) > 0]

    def out_of_bound(self, particle):
        return np.logical_or(np.any(particle.positions < 0, axis=-1),
                             np.any(particle.positions > self.spat_mesh.size, axis=-1))

    def generate_new_particles(self):
        for src in self.particle_sources:
            particles = src.generate_each_step()
            if len(particles.ids):
                particles.ids = self.generate_particle_ids(len(particles.ids))
                self.particle_arrays.append(particles)
        self.shift_new_particles_velocities_half_time_step_back()

    def generate_particle_ids(self, num_of_particles):
        range_of_ids = range(self.max_id + 1, self.max_id + num_of_particles + 1)
        self.max_id += num_of_particles
        return np.array(range_of_ids)

    #
    # Update time grid
    #

    def update_time_grid(self):
        self.time_grid.update_to_next_step()

    #
    # Write domain to file
    #

    def write_step_to_save(self):
        current_step = self.time_grid.current_node
        step_to_save = self.time_grid.node_to_save
        if (current_step % step_to_save) == 0:
            self.write()

    def write(self):
        file_name_to_write = self.construct_output_filename(
            self._output_filename_prefix, self.time_grid.current_node,
            self._output_filename_suffix)
        h5file = h5py.File(file_name_to_write, mode="w")
        if not h5file:
            print("Error: can't open file " + file_name_to_write + \
                  "to save results of initial field calculation!")
            print("Recheck \'output_filename_prefix\' key in config file.")
            print("Make sure the directory you want to save to exists.")
            print("Writing initial fields to file " + file_name_to_write)
        print("Writing step {} to file {}".format(self.time_grid.current_node, file_name_to_write))
        self.save_h5(h5file)
        h5file.close()

    @staticmethod
    def construct_output_filename(output_filename_prefix,
                                  current_time_step,
                                  output_filename_suffix):
        filename = output_filename_prefix + \
                   "{:07d}".format(current_time_step) + \
                   output_filename_suffix
        return filename

    #
    # Free domain
    #

    def free(self):
        print("TODO: free domain.\n")

    #
    # Various functions
    #

    def eval_and_write_fields_without_particles(self):
        self.spat_mesh.clear_old_density_values()
        self.eval_potential_and_fields()
        file_name_to_write = self._output_filename_prefix + "fieldsWithoutParticles" + self._output_filename_suffix
        h5file = h5py.File(file_name_to_write, mode="w")
        if not h5file:
            print("Error: can't open file " + file_name_to_write + \
                  "to save results of initial field calculation!")
            print("Recheck \'output_filename_prefix\' key in config file.")
            print("Make sure the directory you want to save to exists.")
            print("Writing initial fields to file " + file_name_to_write)
        h5file.attrs['class'] = self.__class__.__name__
        self._save_value(h5file, "spat_mesh", self.spat_mesh)
        self._save_value(h5file, "electric_fields", self.electric_fields)
        self._save_value(h5file, "magnetic_fields", self.magnetic_fields)
        self._save_value(h5file, "inner_regions", self.inner_regions)
        h5file.close()

max_stars_repo_name: tnakaicode/ChargedPaticle-LowEnergy

max_issues_repo_name: tnakaicode/ChargedPaticle-LowEnergy

max_forks_repo_name: tnakaicode/ChargedPaticle-LowEnergy