text:
from typing import Callable, Iterable

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation
from matplotlib.axes import Axes
from matplotlib.backend_bases import MouseEvent
from matplotlib.figure import Figure
from matplotlib.lines import Line2D
from matplotlib.widgets import Slider
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from scipy.integrate import odeint

np.warnings.filterwarnings('ignore')

class _MPLPlate(object):
    '''
    This class is a background plate used for all children
    class handlers.

    Parameters
    ----------
        fig : A matplotlib.figure.Figure instance.
        ax : A matplotlib.axes.Axes instance.
    '''

    def __init__(self, fig:Figure, ax:Axes):
        self.fig = fig
        self.ax = ax
        self.ax.set_aspect('equal')

        self.xlim = self.ax.get_xlim()
        self.ylim = self.ax.get_ylim()

        bbox = {'boxstyle': 'round', 'fc': 'w', 'pad': 0.4, 'alpha': 0.7}
        self.annotation = self.ax.annotate('', xy=(0, 0), xytext=(10, 10), textcoords='offset points', bbox=bbox)
        self.annotation.set_visible(False)

class _VectorEventHandler(_MPLPlate):

    def __init__(self, fig:Figure, ax:Axes, x0:float, y0:float, u:float, v:float, name:str, trace_scalars:bool, interactive:bool, **quiver_kwargs):

        super().__init__(fig, ax)

        self.x0 = x0
        self.y0 = y0
        self.x1 = x0 + u
        self.y1 = y0 + v
        self.u = u
        self.v = v
        self.quiver_kwargs = quiver_kwargs

        self.quiver = self.ax.quiver(self.x0, self.y0, self.u, self.v, **self.quiver_kwargs)

        self.trace_state = trace_scalars

        if self.trace_state:
            self.xtrace = self.ax.plot((self.x0, self.x1), (self.y0, self.y0), linestyle='--', color='C0')
            self.ytrace = self.ax.plot((self.x1, self.x1), (self.y0, self.y1), linestyle='--', color='C1')

        if interactive:
            self.dragging_point = False
            self.dragging_vector = False

            self.drag_point = self.ax.scatter(self.x1, self.y1, color='grey', alpha=0.7)
            self.annotation_text = '$\\angle x = {:.2f}^c, \\angle y = {:.2f}^c$\n$x_0 = {:.2f}, y_0 = {:.2f}$\n$x_1 = {:.2f}, y_1 = {:.2f}$\n$u = {:.2f}, v = {:.2f}$\n$mag(\\vec{{' + name + '}}) = {:.2f}$'

            self.fig.canvas.mpl_connect('button_press_event', self.on_click)
            self.fig.canvas.mpl_connect('button_release_event', self.on_release)
            self.fig.canvas.mpl_connect('motion_notify_event', self.on_motion)

    def update_quiver(self, x0:float, y0:float, u:float, v:float) -> None:
        '''
        Updates the drawn quiver and subsequently re-draws the figure.

        Parameters
        ----------
            x0 : x-coordinate of the vector.
            y0 : y-coordinate of the vector.
            u : u scalar of the vector.
            v : v scalar of the vector.

        Notes
        -----
            The quiver is removed from the figure entirely and re-drawn. The
            .set_UVC() method does not work here, since the vector's position
            will be changing, and .set_UVC() only configures the u and v
            scalars.
        '''
        self.drag_point.set_offsets([x0 + u, y0 + v])
        self.quiver.remove()
        self.quiver = self.ax.quiver(x0, y0, u, v, **self.quiver_kwargs)
        self.fig.canvas.draw_idle()

    def update_annotation(self, xdata:float, ydata:float) -> None:
        '''
        Updates the annotation with the given x and y data from the mouse pointer's location.

        Parameters
        ----------
            xdata : Mouse pointer's x position.
            ydata : Mouse pointer's y position.
        '''
        self.annotation.xy = (xdata, ydata)
        self.annotation.set_text(
            self.annotation_text.format(
                np.arctan(self.v / self.u),
                np.arctan(self.u / self.v),
                self.x0,
                self.y0,
                self.x1,
                self.y1,
                self.u,
                self.v,
                np.linalg.norm([self.u, self.v]),
            )
        )
        self.annotation.set_visible(True)
        self.fig.canvas.draw_idle()

    def on_click(self, event:MouseEvent) -> None:
        '''
        Calculates new scalars and initial positions in the event of a drag.

        Parameters
        ----------
            event : Mouse event object emitted by ~Figure.canvas.mpl_connect().
        '''
        if self.drag_point.contains(event)[0] and event.inaxes == self.ax:
            if event.button == 1:
                self.dragging_point = True
            elif event.button == 3:
                self.dragging_vector = True
            self.update_annotation(event.xdata, event.ydata)

    def on_motion(self, event:MouseEvent) -> None:
        '''
        Calculates new scalars and initial positions in the event of a drag.

        Parameters
        ----------
          event : Mouse event object emitted by ~Figure.canvas.mpl_connect().
        '''
        if event.inaxes == self.ax and event.button != 2 and (self.dragging_point or self.dragging_vector):
            self.x1 = event.xdata
            self.y1 = event.ydata

            if self.dragging_point:
                self.u = self.x1 - self.x0
                self.v = self.y1 - self.y0
            elif self.dragging_vector:
                self.x0 = self.x1 - self.u
                self.y0 = self.y1 - self.v

            if self.trace_state:
                self.xtrace[0].set_data((self.x0, self.x1), (self.y0, self.y0))
                self.ytrace[0].set_data((self.x1, self.x1), (self.y0, self.y1))

            self.update_annotation(event.xdata, event.ydata)
            self.update_quiver(self.x0, self.y0, self.u, self.v)

    def on_release(self, event:MouseEvent) -> None:
        '''
        Halts all warping/shifting of the vector if the mouse button has been released. Additionally
        hides annotation.

        Parameters
        ----------
            event : Mouse event object emitted by ~Figure.canvas.mpl_connect().
        '''
        if event.inaxes == self.ax:
            self.annotation.set_visible(False)
            self.fig.canvas.draw_idle()

            if self.drag_point.contains(event)[0] and event.button != 2:
                self.dragging_point = False
                self.dragging_vector = False

                self.x1 = event.xdata
                self.y1 = event.ydata
                self.u = self.x1 - self.x0
                self.v = self.y1 - self.y0

                self.update_quiver(self.x0, self.y0, self.u, self.v)

class _VectorFieldEventHandler(_MPLPlate):

    def __init__(self, fig:Figure, ax:Axes, u:Callable[[float, float], float], v:Callable[[float, float], float],
                 mag:Callable[[float, float], float], div:Callable[[float, float], float], curl:Callable[[float, float], float],
                 name:str, mscale:float, density:int, normalize:bool, colorbar:bool, interactive:bool, **quiver_kwargs):

        super().__init__(fig, ax)

        self.u = u
        self.v = v

        self.mag = mag
        self.div = div
        self.curl = curl

        self.scale = mscale
        self.quiver_kwargs = quiver_kwargs

        self.normalize_state = normalize
        self.colorbar_state = colorbar

        components = self.create_field(self.scale, density)

        self.quiver = self.ax.quiver(*components, **self.quiver_kwargs)

        if self.colorbar_state:
            plt.subplots_adjust(left=0.05)
            c = components[-1]

            if self.normalize_state:
                self.colorbar_label = 'Sampled Magnitude (Scale: {:.2f}, Normalized)'
            else:
                self.colorbar_label = 'Sampled Magnitude (Scale: {:.2f})'

            cax = inset_axes(
                self.ax,
                width='5%',
                height='100%',
                loc='lower left',
                bbox_to_anchor=(1.05, 0., 1, 1),
                bbox_transform=self.ax.transAxes,
                borderpad=0,
                axes_kwargs={'zorder': -1}
            )

            self.colorbar = plt.colorbar(
                self.quiver,
                cax=cax,
                cmap=self.quiver_kwargs['cmap'],
                label=self.colorbar_label.format(self.scale),
                ticks=np.linspace(c.min(), c.max(), 5)
            )

        if interactive:
            self.annotation_text = '$x = {:.2f}, y = {:.2f}$\n$mag(\\vec{{' + name + '}})(x, y) = {:.2f}$\n$div(\\vec{{' + name + '}})(x, y) = {:.2f}$\n$curl(\\vec{{' + name + '}})(x, y) = {:.2f}$'

            self.fig.canvas.mpl_connect('button_press_event', self.on_click)
            self.fig.canvas.mpl_connect('button_release_event', self.on_release)

            scale_range = round(max(abs(val) for val in (self.xlim + self.ylim)) / 4)
            if scale_range < 1: scale_range = 1

            divider = make_axes_locatable(self.ax)

            scale_ax = divider.append_axes('bottom', size='3%', pad=0.6)
            self.scale_slider = Slider(scale_ax, 'Scale', -scale_range, scale_range, valinit=mscale)
            self.scale_slider.on_changed(self.slider_update)

            density_ax = divider.append_axes('bottom', size='3%', pad=0.1)
            self.density_slider = Slider(density_ax, 'Density', 4, 100, valinit=density, valstep=1)
            self.density_slider.on_changed(self.slider_update)

    def on_release(self, event:MouseEvent) -> None:
        '''
        Clears annotation on release of the mouse button.

        Parameters
        ----------
            event : Mouse event emitted by ~Figure.canvas.mpl_connect()
        '''
        if event.inaxes == self.ax:
            self.annotation.set_visible(False)
            self.fig.canvas.draw_idle()

    def on_click(self, event:MouseEvent) -> None:
        '''
        Displays and updates annotation on mouse click.

        Parameters
        ----------
            event : Mouse event emitted by ~Figure.canvas.mpl_connect()
        '''
        if event.inaxes == self.ax:
            self.annotation.xy = (event.xdata, event.ydata)
            self.annotation.set_text(
                self.annotation_text.format(
                    event.xdata,
                    event.ydata,
                    self.mag(event.xdata, event.ydata),
                    self.div(event.xdata, event.ydata),
                    self.curl(event.xdata, event.ydata),
                    )
                )
            self.annotation.set_visible(True)
            self.fig.canvas.draw_idle()

    def slider_update(self, _) -> None:
        '''
        Updates and re-draws field plot based on the scale and density slider values.
        '''

        components = self.create_field(self.scale_slider.val, self.density_slider.val)

        self.quiver.remove()
        self.quiver = self.ax.quiver(*components, **self.quiver_kwargs)

        if self.colorbar_state:
            self.colorbar.set_label(self.colorbar_label.format(self.scale_slider.val))

        self.fig.canvas.draw_idle()

    def create_field(self, scale:float, density:int) -> tuple:
        '''
        Calculates all vector and vector positions with a given scale and density.

        Parameters
        ----------
            scale : Scale value applied to vector lengths.
            density : Density of the field.

        Returns
        -------
            tuple : A tuple of ndarrays, consisting of the x, y, u, v and c
                    values of the field respectively.
        '''

        x, y = np.meshgrid(
            np.linspace(*self.xlim, int(density)),
            np.linspace(*self.ylim, int(density))
        )

        u, v = (
            self.u(x, y),
            self.v(x, y)
        )

        m = np.sqrt(u**2 + v**2)

        if self.normalize_state:
            with np.errstate(all='ignore'):
                u = (u / m) * scale
                v = (v / m) * scale
        else:
            u = u * scale
            v = v * scale

        return x, y, u, v, m

class _ParticleSimulationHandler(_MPLPlate):

    def __init__(self, fig:Figure, ax:Axes, u:Callable[[float, float], float], v:Callable[[float, float], float], pts:Iterable[tuple],
                 frames:int, dt:float, blit:bool, fmt:str, color:str, alpha:float, **plot_kwargs):

        super().__init__(fig, ax)

        self.u = u
        self.v = v

        self.ln, = self.ax.plot([], [], fmt, color=color, alpha=alpha, **plot_kwargs)

        if pts is None:
            pts = np.array((np.random.uniform(*self.xlim, 50), np.random.uniform(*self.ylim, 50))).transpose()

        self.ani = FuncAnimation(self.fig, self.particle_update, interval=1, frames=frames, blit=blit, fargs=(pts, dt))

    def solve_ode(self, f:Callable, pts:Iterable[tuple], dt:float) -> list:
        '''
        Solves for the displacement of each particle with respect to the change in time.

        Parameters
        ----------
            f : Function integrand. In our case, the u and v scalar functions.
            pts : Array of coordinate pairs.
            dt : The change in time from one frame to the next.

        Returns
        -------
            list : A list of all updated particle positions.
        '''
        return [odeint(f, pt, [0, dt])[-1] for pt in pts]

    def get_vels(self, pt:tuple, _) -> list:
        '''
        Calculates the velocity of the particle at a specific point.

        Parameters
        ----------
            pt : An x, y coordinate pair.
            _ : Dummy time parameter required for odeint().
                - Unused since scalar functions don't depend on time.

        Returns
        -------
            list : The x and y velocity at the given point (index 0 and index 1, respectively).
        '''
        return [self.u(*pt), self.v(*pt)]

    def remove_pts(self, pts:np.ndarray) -> np.ndarray:
        '''
        Removes points that are outside the bounds of the axes.

        Parameters
        ----------
            pts : Array of coordinate pairs.

        Returns
        -------
            numpy.ndarray : The same array but with out-of-boundary points removed.

        Credits
        -------
            Tony S. Yu, Ph.D.
        '''
        if len(pts) == 0: return []

        out_x = (pts[:, 0] < self.xlim[0]) | (pts[:, 0] > self.xlim[1])
        out_y = (pts[:, 1] < self.ylim[0]) | (pts[:, 1] > self.ylim[1])
        keep = ~(out_x | out_y)

        return pts[keep]

    def particle_update(self, frame:int, *fargs) -> Line2D:
        '''
        Calculates particle displacements, removes out-of-axes points, and subsequently updates the
        axes with new particle positions.

        Parameters
        ----------
            frame : Current frame of the animation.
            *fargs : A list of the intial particle positions at index 0 (type: Iterable), and the dt value at index 1 (type: float).

        Returns
        -------
            matplotlib.lines.Line2D : The updated Line2D array.
        '''
        # self.pts needs to be a data member, though it may not look like it at first.
        # transition arrays would be reset back to the initial if not a data member.

        if frame == 0: self.pts = fargs[0]

        self.pts = np.asarray(self.solve_ode(self.get_vels, self.pts, fargs[1]))
        self.pts = np.asarray(self.remove_pts(self.pts))

        self.pts.shape = (self.pts.shape[0], 2)  # rudimentary way of ensuring .transpose() runs properly
        x, y = self.pts.transpose()

        self.ln.set_data(x, y)

        return self.ln,

max_stars_repo_name: braedynl/VectorGraphingToolkit

max_issues_repo_name: braedynl/VectorGraphingToolkit

max_forks_repo_name: braedynl/VectorGraphingToolkit


text:
import itertools
import numpy as np
from typing import Any, Dict, Iterable, Iterator, List, Optional, Set, Tuple, Union

from pandemicsim.agent import Agent

GridContent = Union[Optional[Agent], Set[Agent]]
FloatCoordinate = Union[Tuple[float, float], np.ndarray]

class ContinuousSpace:
	""" Continuous space where each agent can have an arbitrary position.

	Assumes that all agents are point objects, and have a pos property storing
	their position as an (x, y) tuple. This class uses a numpy array internally
	to store agent objects, to speed up neighborhood lookups.

	"""

	# _grid = None

	def __init__(self, width: float, height: float) -> None:
		""" Create a new continuous space."""

		self.width = width
		self.height = height

		self._agent_points = None	# type: np.ndarray of FloatCoordinates
		self._index_to_agent = {}  # type: Dict[int, Agent]
		self._agent_to_index = {}  # type: Dict[Agent, int]

	def place_agent(self, agent: Agent, pos: FloatCoordinate) -> None:
		""" Place a new agent in the space.

		Args:
			agent: Agent object to place.
			pos: Coordinate tuple for where to place the agent.

		"""
		if self._agent_points is None:
			self._agent_points = np.array([pos])
		else:
			self._agent_points = np.append(self._agent_points, np.array([pos]), axis=0)
		self._index_to_agent[self._agent_points.shape[0] - 1] = agent
		self._agent_to_index[agent] = self._agent_points.shape[0] - 1
		agent.pos = pos

	def move_agent(self, agent: Agent, pos: FloatCoordinate) -> None:
		""" Move an agent from its current position to a new position.

		Args:
			agent: The agent object to move.
			pos: Coordinate tuple to move the agent to.

		"""
		idx = self._agent_to_index[agent]
		self._agent_points[idx, 0] = pos[0]
		self._agent_points[idx, 1] = pos[1]
		agent.pos = pos

	def remove_agent(self, agent: Agent) -> None:
		""" Remove an agent from the simulation.

		Args:
			agent: The agent object to remove
			"""
		if agent not in self._agent_to_index:
			raise Exception("Agent does not exist in the space")
		idx = self._agent_to_index[agent]
		del self._agent_to_index[agent]
		max_idx = max(self._index_to_agent.keys())
		# Delete the agent's position and decrement the index/agent mapping
		self._agent_points = np.delete(self._agent_points, idx, axis=0)
		for a, index in self._agent_to_index.items():
			if index > idx:
				self._agent_to_index[a] = index - 1
				self._index_to_agent[index - 1] = a
		# The largest index is now redundant
		del self._index_to_agent[max_idx]
		agent.pos = None

	def get_neighbors(
		self, pos: FloatCoordinate, radius: float, include_center: bool = True
	) -> List[GridContent]:
		""" Get all objects within a certain radius.

		Args:
			pos: (x,y) coordinate tuple to center the search at.
			radius: Get all the objects within this distance of the center.
			include_center: If True, include an object at the *exact* provided
							coordinates. i.e. if you are searching for the
							neighbors of a given agent, True will include that
							agent in the results.

		"""
		deltas = np.abs(self._agent_points - np.array(pos))
		dists = deltas[:, 0] ** 2 + deltas[:, 1] ** 2

		(idxs,) = np.where(dists <= radius ** 2)
		neighbors = [
			self._index_to_agent[x] for x in idxs if include_center or dists[x] > 0
		]
		return neighbors

	def get_distance(self, pos_1: FloatCoordinate, pos_2: FloatCoordinate) -> float:
		""" Get the distance between two point, accounting for toroidal space.

		Args:
			pos_1, pos_2: Coordinate tuples for both points.

		"""
		x1, y1 = pos_1
		x2, y2 = pos_2

		dx = np.abs(x1 - x2)
		dy = np.abs(y1 - y2)

		return np.sqrt(dx * dx + dy * dy)

	def out_of_bounds(self, pos: FloatCoordinate) -> bool:
		""" Check if a point is out of bounds. """
		x, y = pos
		return x < self.x_min or x >= self.x_max or y < self.y_min or y >= self.y_max

	# def plot(self):
	# 	for point in self._agent_points:
	# 		plt.plot(point[0], point[1])

ext: py

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Feb  4 19:26:53 2019

@author: MAGESHWARAN
"""

import numpy as np
import math
from scipy import optimize as op

class Regularization(object):

    def __init__(self, alpha):
        self.alpha = alpha

    def __call__(self, weights):
        return 0.5 * self.alpha * weights.T.dot(weights)

    def grad(self, weights):
        return self.alpha * weights.T.dot(weights)

class LogRegWithRegularization(object):

    def __init__(self, n_iterations = 10000, use_optimizer = True, reg_param = 0.01, learning_rate = 0.001):
        self.n_iterations = n_iterations
        self.learning_rate = learning_rate
        self.regularization = Regularization(alpha = reg_param)
        self.use_optimizer = use_optimizer

    def initialize_weights(self, n_features):
        # initialize the weights based on number of features
        limit = 1 / math.sqrt(n_features)
        weights = np.random.uniform(-limit, limit, (n_features,))
        return weights

    # using sigmoid function for prediction
    def sigmoid(self, X, weights):
        z = X.dot(weights)
        h = 1 / (1 + np.exp(-z))
        return h

    def loss(self, weights, X, y):
       y_pred = self.sigmoid(X, weights)
       # calculating the cost / loss of the model
       J = np.mean(-(y * np.log(y_pred)) - ((1 - y) * np.log(1 - y_pred)) + self.regularization(self.weights))
       return J

    def gradient_descent(self, weights, X, y):
        y_pred = self.sigmoid(X, weights)
        # Calculate gradient
        grad = - (y - y_pred).dot(X) + self.regularization.grad(weights)
        # return grad / X.shape[0]
        return grad

    # Train the model
    def fit(self, X, y):
        # add bias units
        X = np.insert(X, 0, 1, axis = 1)

        self.training_errors = []
        self.weights = self.initialize_weights(n_features = X.shape[1])

        # Gradient Descent using n_iterations
        if not self.use_optimizer:
            for i in range(self.n_iterations):
                # Find loss / training error
                mse = self.loss(self.weights, X, y)
                self.training_errors.append(mse)
                grad = self.gradient_descent(self.weights, X, y)
                # gradient update
                self.weights -= self.learning_rate * grad
        else:
            result = op.minimize(fun = self.loss, x0 = self.weights, args = (X, y),
                                       method = "TNC", jac = self.gradient_descent)
            self.weights = result.x

    def predict(self, X):
        X = np.insert(X, 0, 1, axis = 1)
        y_pred = self.sigmoid(X, self.weights)
        y_pred = np.round(y_pred).astype(int)
        return y_pred

max_stars_repo_name: MageshDominator/data-science-py

max_issues_repo_name: MageshDominator/data-science-py

max_forks_repo_name: MageshDominator/data-science-py

reason: import numpy,from scipy


text:
'''
By adidinchuk. adidinchuk@gmail.com.
https://github.com/adidinchuk/tf-neural-net
'''

import hyperparams as hp
import data as d
import numpy as np
import network as nwk

import os
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'

data = d.import_data('promoted.csv', headers=True)

features_numeric = d.zero_to_mean([[float(row[2].strip()) if row[2].strip() else 0.,
                     float(row[3].strip()) if row[3].strip() else 0.,
                     float(row[5].strip()) if row[5].strip() else 0.]
                    for row in data])
# normalize numeric features
features_numeric = np.transpose(d.normalize(features_numeric))

features_string = [[row[6].strip() if row[6].strip() else 'N/A',
                    row[7].strip() if row[7].strip() else 'N/A']
                   for row in data]

# expand the categorical columns into binary categories
row4_expansion, row4_translation = d.expand_categorical_feature(np.transpose(features_string)[0])
row4_expansion = row4_expansion
print('Row 4 expanded into the following categories: ' + str(row4_translation))
row5_expansion, row5_translation = d.expand_categorical_feature(np.transpose(features_string)[1])
row5_expansion = row5_expansion
print('Row 5 expanded into the following categories: ' + str(row5_translation))

# concatenate numerical data with categorical features
features = np.append(features_numeric, row4_expansion, axis=0)
features = np.append(features, row5_expansion, axis=0)
features = np.transpose(features)

# extract target data
target = np.array([[float(row[1])] for row in data])

# split into test and train data sets
train_size = round(len(features) * (1 - hp.validation_size))
train_index = np.random.choice(len(features), size=train_size)
test_index = list(set(range(len(features))) - set(train_index))

# training data
train_inputs = features[train_index]
train_outputs = target[train_index]

# testing data
test_inputs = features[test_index]
test_outputs = target[test_index]

nn = nwk.Network(len(features[0]), len(target[0]), hp.nn_structure, hp.nn_activations)
nn.train(train_inputs, train_outputs, test_inputs, test_outputs, auto_balance=hp.auto_balance, plot=True,
         loss_function=hp.loss_function, batch_size=hp.batch_size, epochs=hp.epochs, learning_rate=hp.learning_rate)

max_stars_repo_name: adidinchuk/tf-neural-net

max_issues_repo_name: adidinchuk/tf-neural-net

max_forks_repo_name: adidinchuk/tf-neural-net


text:
import json
import os
import warnings

import casadi as ca
import numpy as np
import pandas as pd
import pandas.testing as pdt
import pytest
from scipy.signal import chirp

from skmid.integrator import RungeKutta4
from skmid.models import DynamicModel
from skmid.models import generate_model_attributes

@pytest.fixture
def load_non_linear_model_data():
    """Generate input signal"""
    CWD = os.getcwd()
    DATA_DIR = "data"
    SUB_DATA_DIR = "non_linear_model"

    U = pd.read_csv(
        filepath_or_buffer=os.path.join(CWD, DATA_DIR, SUB_DATA_DIR, "u_data.csv"),
        index_col=0,
    )
    Y = pd.read_csv(
        filepath_or_buffer=os.path.join(CWD, DATA_DIR, SUB_DATA_DIR, "y_data.csv"),
        index_col=0,
    )

    # reading the data from the file
    with open(
        os.path.join(CWD, DATA_DIR, SUB_DATA_DIR, "settings.json"), mode="r"
    ) as j_object:
        settings = json.load(j_object)

    return (U, Y, settings)

# test f(x,u) with multiple input, multiple size, list
@pytest.fixture
def generate_input_signal():
    """Generate input signal"""

    N = 2000  # Number of samples
    fs = 500  # Sampling frequency [hz]
    t = np.linspace(0, (N - 1) * (1 / fs), N)
    df_input = pd.DataFrame(
        data={
            "chirp": 2 * chirp(t, f0=1, f1=10, t1=5, method="logarithmic"),
            "noise": 2 * np.random.random(N),
        },
        index=t,
    )

    return (df_input, fs)

@pytest.fixture
def generate_step_signal():
    """Generate Step input signal"""

    N = 1000  # Number of samples
    fs = 25  # Sampling frequency [hz]
    t = np.linspace(0, (N - 1) * (1 / fs), N)
    df_input = pd.DataFrame(index=t).assign(
        step=lambda x: np.where(x.index < t[int(N / 4)], 0, 1)
    )

    return (df_input, fs)

@pytest.fixture
def generate_inpulse_signal():
    """Generate Impulse input signal"""

    N = 500  # Number of samples
    fs = 50  # Sampling frequency [hz]
    t = np.linspace(0, (N - 1) * (1 / fs), N)
    df_input = pd.DataFrame(index=t).assign(inpulse=np.zeros(N))
    df_input.iloc[0] = 1

    return (df_input, fs)

class TestRungeKutta4:
    """Test class for function generate_model_parameters."""

    def test_model_with_states(self):
        """Test simulation with model dx=f(x)."""

        (x, _, _) = generate_model_attributes(
            state_size=1, input_size=0, parameter_size=0
        )

        # initialize first-order model
        tau = 1
        sys = DynamicModel(state=x, model_dynamics=[-(1 / tau) * x])

        n_steps = 50
        x0 = [1]
        rk4 = RungeKutta4(model=sys)
        _ = rk4.simulate(initial_condition=x0, n_steps=n_steps)

        df_X = rk4.state_sim_
        df_Y = rk4.output_sim_

        # check equality of dataframe
        pdt.assert_frame_equal(df_X, df_Y)

        # check size of dataframe. Note: +1 is because it includes the initial condition
        assert len(df_X) == n_steps + 1
        assert len(df_Y) == n_steps + 1

        # check frequency is assigned correctly fs=1 default
        assert all(np.diff(df_Y.index) == 1)

        # check no missing values
        assert df_X.notna().all().values
        assert df_Y.notna().all().values

        # check model consistency: convergence with no bias
        assert df_X.iloc[0].values == pytest.approx(x0)
        assert df_X.iloc[-1].values == pytest.approx(0)

    def test_model_with_states_input(self, generate_step_signal):
        """Test simulation with model dx=f(x,u)."""

        (df_input, fs) = generate_step_signal

        # initialize first-order model
        (x, u, _) = generate_model_attributes(
            state_size=1, input_size=1, parameter_size=0
        )

        tau, kp = 1, 1
        sys = DynamicModel(
            state=x, input=u, model_dynamics=[-(1 / tau) * x + (kp / tau) * u]
        )

        rk4 = RungeKutta4(model=sys, fs=fs)
        _ = rk4.simulate(initial_condition=[0], input=df_input)

        df_X = rk4.state_sim_
        df_Y = rk4.output_sim_

        # check equality of dataframe
        pdt.assert_frame_equal(df_X, df_Y)

        # check size of dataframe
        assert len(df_X) == len(df_input) + 1
        assert len(df_Y) == len(df_input) + 1

        # check no missing values
        assert df_X.notna().all().values
        assert df_Y.notna().all().values

        # check model consistency: convergence with no bias
        assert df_input.iloc[-1].values == pytest.approx(df_X.iloc[-1].values)

    def test_model_with_states_parameter(self):
        """Test simulation with model dx=f(x,u)."""

        (x, _, tau) = generate_model_attributes(
            state_size=1, input_size=0, parameter_size=1
        )

        sys = DynamicModel(state=x, parameter=tau, model_dynamics=[-(1 / tau) * x])

        n_steps = 50
        x0 = [1]
        tau_num = [1]
        rk4 = RungeKutta4(model=sys)
        _ = rk4.simulate(initial_condition=x0, parameter=tau_num, n_steps=n_steps)

        df_X = rk4.state_sim_
        df_Y = rk4.output_sim_

        # check equality of dataframe
        pdt.assert_frame_equal(df_X, df_Y)

        # check size of dataframe. Note: +1 is because it includes the initial condition
        assert len(df_X) == n_steps + 1
        assert len(df_Y) == n_steps + 1

        # check frequency is assigned correctly fs=1 default
        assert all(np.diff(df_Y.index) == 1)

        # check no missing values
        assert df_X.notna().all().values
        assert df_Y.notna().all().values

        # check model consistency: convergence to 1 with no bias
        assert df_X.iloc[0].values == pytest.approx(x0)
        assert df_X.iloc[-1].values == pytest.approx(0)

    def test_model_with_states_input_parameter(self, generate_step_signal):
        """Test simulation with model dx=f(x,u,p)."""

        (df_input, fs) = generate_step_signal

        (x, u, parameter) = generate_model_attributes(
            state_size=1, input_size=1, parameter_size=2
        )

        tau, kp = parameter[0], parameter[1]
        sys = DynamicModel(
            state=x,
            input=u,
            parameter=parameter,
            model_dynamics=[-(1 / tau) * x + (kp / tau) * u],
        )

        parameter_num = [1, 1]
        rk4 = RungeKutta4(model=sys, fs=fs)
        _ = rk4.simulate(initial_condition=[0], input=df_input, parameter=parameter_num)

        df_X = rk4.state_sim_
        df_Y = rk4.output_sim_

        # check equality of dataframe
        pdt.assert_frame_equal(df_X, df_Y)

        # check size of dataframe
        assert len(df_X) == len(df_input) + 1
        assert len(df_Y) == len(df_input) + 1

        # check no missing values
        assert df_X.notna().all().values
        assert df_Y.notna().all().values

        # check model consistency: convergence with no bias
        assert df_input.iloc[-1].values == pytest.approx(df_X.iloc[-1].values)

    def test_different_format_input(self, generate_inpulse_signal):
        """."""

        (df_input, fs) = generate_inpulse_signal

        (x, u, _) = generate_model_attributes(
            state_size=1, input_size=1, parameter_size=0
        )

        tau, kp = 1, 1
        sys = DynamicModel(
            state=x, input=u, model_dynamics=[-(1 / tau) * x + (kp / tau) * u]
        )

        rk4 = RungeKutta4(model=sys, fs=fs)

        input_dict = {
            "dataframe": df_input,
            "dataframe_with_non_consistent_fs": df_input.reset_index(drop=True),
            "series": df_input.iloc[:, 0],
            "numpy_array": df_input.values,
            "numpy_array_1dim": df_input.values.squeeze(),
        }

        for format, input_ in input_dict.items():
            message = f"Test failed with input format {format}"

            _ = rk4.simulate(initial_condition=[0], input=input_)

            df_X = rk4.state_sim_
            df_Y = rk4.output_sim_

            # check equality of dataframe
            assert isinstance(df_X, pd.DataFrame), message
            assert isinstance(df_Y, pd.DataFrame), message

            pdt.assert_frame_equal(df_X, df_Y)

            # check size of dataframe
            assert len(df_X) == len(df_input) + 1, message
            assert len(df_Y) == len(df_input) + 1, message

            # check no missing values
            assert df_X.notna().all().values, message
            assert df_Y.notna().all().values, message

            # check model consistency: convergence with no bias
            assert df_input.iloc[-1].values == pytest.approx(
                df_X.iloc[-1].values, rel=1e-7, abs=1e-6
            ), message

    def test_mimo_system(self):

        # generate dirac impulse
        N = 3000  # Number of samples
        fs = 50  # Sampling frequency [hz]
        t = np.linspace(0, (N - 1) * (1 / fs), N)
        df_input = pd.DataFrame(
            data=np.zeros(shape=(N, 2)), columns=["u1", "u2"], index=t
        )
        df_input.iloc[0, :] = [0.1, -0.1]

        (x, u, param) = generate_model_attributes(
            state_size=2, input_size=2, parameter_size=2
        )

        x1, x2 = x[0], x[1]
        u1, u2 = u[0], u[1]
        k1, k2 = param[0], param[1]

        # chemical reactor: mimo
        rhs = [u1 - k1 * x1, u1 * u2 / x1 - u1 * x2 / x1 - k2 * x2]

        sys = DynamicModel(
            state=x,
            input=u,
            parameter=param,
            model_dynamics=rhs,
            parameter_name=["k1", "k2"],
        )

        rk4 = RungeKutta4(model=sys, fs=fs)

        input_dict = {
            "dataframe": df_input,
            "dataframe_with_non_consistent_fs": df_input.reset_index(drop=True),
            "numpy_array": df_input.values,
        }

        for format, input_ in input_dict.items():
            message = f"Test failed with input format {format}"

            _ = rk4.simulate(
                initial_condition=[-0.15, 0.1], input=input_, parameter=[0.1, 0.5]
            )

            df_X = rk4.state_sim_
            df_Y = rk4.output_sim_

            # check equality of dataframe
            assert isinstance(df_X, pd.DataFrame), message
            assert isinstance(df_Y, pd.DataFrame), message

            pdt.assert_frame_equal(df_X, df_Y)

            # check size of dataframe
            assert len(df_X) == len(df_input) + 1, message
            assert len(df_Y) == len(df_input) + 1, message

            # check no missing values
            assert df_X.notna().all().all(), message
            assert df_Y.notna().all().all(), message

            # check model consistency: convergence with no bias
            assert df_input.iloc[-1].values == pytest.approx(
                df_X.iloc[-1].values, rel=1e-3, abs=1e-3
            ), message

    def test_input_sample_frequency(self):

        (x, u, _) = generate_model_attributes(
            state_size=1, input_size=1, parameter_size=0
        )

        tau, kp = 1, 1
        sys = DynamicModel(
            state=x, input=u, model_dynamics=[-(1 / tau) * x + (kp / tau) * u]
        )

        with pytest.raises(ValueError):
            _ = RungeKutta4(model=sys, fs=-12)

            _ = RungeKutta4(model=sys, fs=0)

            _ = RungeKutta4(model=sys, fs=50, n_steps_per_sample=0)

            _ = RungeKutta4(model=sys, fs=50, n_steps_per_sample=0.1)

            _ = RungeKutta4(model=sys, fs=50, n_steps_per_sample=-np.sqrt(0.1))

max_stars_repo_name: G-Licitra/scikit-sysid

max_issues_repo_name: G-Licitra/scikit-sysid

max_forks_repo_name: G-Licitra/scikit-sysid


text:
import numpy as np
import os
import pickle

def get_baseline_estimates(ratings):
    #print('inside get_baseline_estimates')
    if not os.path.isfile("baseline_estimates.pkl"):
        user_count = ratings.shape[0]
        movie_count = ratings.shape[1]

        avg_user = []
        avg_movie = []

        global_mean = 0

        ratings_count = 0
        for x, row in enumerate(ratings):
            nonzero_ratings = 0
            sum_ratings = 0
            for rating in row:
                if rating != 0:
                    nonzero_ratings += 1
                    sum_ratings += rating

            global_mean += sum_ratings
            ratings_count += nonzero_ratings
            if nonzero_ratings == 0:
                avg_user.append(0)
            else:
                avg_user.append(float(sum_ratings) / nonzero_ratings)

        global_mean /= ratings_count

        for i, column in enumerate(np.transpose(ratings)):
            nonzero_ratings = 0
            sum_ratings = 0
            for rating in column:
                if rating != 0:
                    nonzero_ratings += 1
                    sum_ratings += rating
            if nonzero_ratings == 0:
                avg_movie.append(0)
            else:
                avg_movie.append(float(sum_ratings) / nonzero_ratings)

        baseline_est = np.zeros(ratings.shape)

        for i in range(user_count):
            for x in range(movie_count):
                if ratings[i, x] != 0:
                    baseline_est[i, x] = avg_movie[x] + avg_user[i] - global_mean

        pickle.dump(baseline_est, open("baseline_estimates.pkl", "wb"))
    else:
        with open('baseline_estimates.pkl', 'rb') as pickle_file:
            baseline_est = pickle.load(pickle_file)

    return baseline_est


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

"""

This is the interface of the MILES input files as follows:

  observations_fitMIR.h5   -- data + new initial parameters
  input_fitMIR_master.h5   -- general param
  input_fitMIR_model.h5    -- model param
  input_fitMIR_extra.h5    -- extra param
  input_fitMIR_analysis.h5 -- param for analysis

"""

import os
import math
import numpy as np

## astylo
from astylo.arrlib import closest
from astylo.iolib import read_fits, write_hdf5, read_hdf5

## local
from utilities import res, TABLine, TABand, partuning

## Path
##------
mroot = os.path.dirname(os.path.abspath(__file__))+'/'
h5_obs = mroot+'dat/observations_fitMIR'
h5_master = mroot+'input_fitMIR_master'
h5_model = mroot+'input_fitMIR_model'
h5_extra = mroot+'input_fitMIR_extra'
h5_analysis = mroot+'input_fitMIR_analysis'

##-----------------------------
## Set parameters
##-----------------------------
noisy = False # verbose/debug for this routine
chi2init = True # True if using Chi2 results as HB init param

##-----------------------------
## Write observations_fitMIR.h5 (part 1)
##-----------------------------
z = 0.0017
fits_obs = mroot+'dat/M83s' # obs
fits_unc = fits_obs+'_unc' # unc
wvl_inf = 1. # min wvl
wvl_sup = 21. # max wvl
x_inf = 58
x_sup = 60
y_inf = 23
y_sup = 25
# spec_unit = 'MKS' # W.m-2.Hz-1.sr-1
spec_unit = 'MJyovsr' # MJy.sr-1

## Read FITS
dset = read_fits(fits_obs, fits_unc)

## Truncate wavelength range
ind_inf = closest(dset.wave, wvl_inf)
ind_sup = closest(dset.wave, wvl_sup)
wave = dset.wave[ind_inf:ind_sup]
wave = wave / (1+z)

data = dset.data[ind_inf:ind_sup,y_inf:y_sup,x_inf:x_sup]
unc = dset.unc[ind_inf:ind_sup,y_inf:y_sup,x_inf:x_sup]
## convert MJy/sr to W/m2/Hz/sr (MKS)
if spec_unit=='MKS':
    data = data * 1.e-20
    unc = unc * 1.e-20

## Mask NaNs
mask = np.isnan(data) * 1

## Write HDF5
##------------
if (not chi2init):
    write_hdf5(h5_obs, 'Spectral unit', [spec_unit], verbose=noisy)
    write_hdf5(h5_obs, 'Wavelength (microns)', wave, append=True, verbose=noisy)
    write_hdf5(h5_obs, 'FnuOBS ('+spec_unit+')', data, append=True, verbose=noisy)
    write_hdf5(h5_obs, 'dFnuOBS ('+spec_unit+')', unc, append=True, verbose=noisy)
    write_hdf5(h5_obs, 'NaN mask', mask, append=True, verbose=noisy)
    write_hdf5(h5_obs, 'Chi2init', ['F'], append=True, verbose=noisy)
else:
    write_hdf5(h5_obs, 'Chi2init', ['T'], append=True, verbose=noisy)

##--------------------------------
## Write input_fitMIR_master.h5
##--------------------------------
verbose = 'T'
Nmcmc = 50
NiniMC = 0
calib = 'F'
robust_RMS = 'F'
robust_cal = 'F'
skew_RMS = 'F'
newseed = 'F'
dostop = 'F'
newinit = 'F'

## Chi2 results are used as HB init param
if (chi2init):
    newinit = 'T'

## Write HDF5
##------------
write_hdf5(h5_master, 'Spectral unit', [spec_unit], verbose=noisy)
write_hdf5(h5_master, 'verbose', [verbose], append=True, verbose=noisy)
write_hdf5(h5_master, 'Nmcmc', [Nmcmc], append=True, verbose=noisy)
write_hdf5(h5_master, 'NiniMC', [NiniMC], append=True, verbose=noisy)
write_hdf5(h5_master, 'calib', [calib], append=True, verbose=noisy)
write_hdf5(h5_master, 'robust_RMS', [robust_RMS], append=True, verbose=noisy)
write_hdf5(h5_master, 'robust_cal', [robust_cal], append=True, verbose=noisy)
write_hdf5(h5_master, 'skew_RMS', [skew_RMS], append=True, verbose=noisy)
write_hdf5(h5_master, 'newseed', [newseed], append=True, verbose=noisy)
write_hdf5(h5_master, 'newinit', [newinit], append=True, verbose=noisy)
write_hdf5(h5_master, 'dostop', [dostop], append=True, verbose=noisy)
# write_hdf5(h5_master, '', , append=True, verbose=noisy)

##--------------------------------
## Write input_fitMIR_model.h5
##--------------------------------
labQ = ['ACH2_Z96             ', # 9
        'BE_Z96               ', # 10
        'Sil_D03              '] # 23
labL = ['H2S7  ', # 3
        'H2S5  ', # 9
        'ArII  ', # 11
        'ArIII1', # 20
        'H2S3  ', # 23
        'HeII2 ', # 24
        'SIV   ', # 26
        'H2S2  ', # 27
        'NeII  ', # 29
        'NeIII1', # 33
        'H2S1  ', # 34
        'SIII1 '] # 36
labB = ['Main 3.3     ', # 1
        'Main 3.4     ', # 2
        'Main 6.2 (1) ', # 7*
        'Main 6.2 (2) ', # 8
        'Plateau 7.7  ', # 12
        'Main 7.7 (1) ', # 13*
        'Main 7.7 (2) ', # 14
        'Main 8.6     ', # 16*
        'Small 11.0   ', # 19
        'Main 11.2    ', # 20*
        'Plateau 11.3 ', # 21
        'Main 12.7 (1)', # 24
        'Main 12.7 (2)', # 25
        'Plateau 17.0 '] # 30

ALline = False
ALband = True

if (ALline):
    labL = []
    for tabL in TABLine:
        if (tabL['wave']>wave[0] and tabL['wave']<wave[-1]):
            labL.append(tabL['label'])
if (ALband):
    labB = []
    for tabB in TABand:
        if (tabB['wave']>wave[0] and tabB['wave']<wave[-1]):
            labB.append(tabB['label'])

## dictune, as input of the partuning function,
## is a list of dict containing param tuning info
dictune = [ dict([ ('name','default'),
                   ('namall','default'),
                   ('fixed','T'),
                   ('limits',(0.,0.)),
                   ('limited',('F','F')),
                   ('model','T'),
                   ('hyper','F'),
                   ('tied',''),
                   ('value',0.) ]),

            ## Add your config below
            ##=======================

            ## Extensive param:
            ## lnMovd2, lnIline, lnIband, lnFstar,
            ##-----------------
            dict([ ('namall','lnMovd2'),
                   ('fixed','F'),
                   # ('limits',(-25.,0.)),
                   # ('limited',('T','F')),
                   # ('hyper','T'),
            ]),

            dict([ ('namall','lnIline'),
                   ('fixed','F'),
                   # ('limits',(-25.,10.)),
                   # ('limited',('T','F')),
                   # ('hyper','T'),
            ]),

            dict([ ('namall','lnIband'),
                   ('fixed','F'),
                   # ('limits',(-25.,10.)),
                   # ('limited',('T','F')),
                   # ('hyper','T'),
            ]),

            dict([ ('name','lnFstar'),
                   ('fixed','F'),
                   # ('limits',(-25.,10.)),
                   # ('limited',('T','F')),
                   # ('hyper','T'),
                   ]),

            ## Intensive param:
            ## lnT, Cline, Cband,
            ## (fixed) Wline, WSband, WLband, lnAv,
            ##-----------------
            dict([ ('namall','lnT'),
                   ('fixed','F'),
                   ('limited',('T','T')),
                   ('limits',(np.log(50.),np.log(500.))), # LOG( (50,500) K )
                   # ('hyper','T'),
            ]),

            # dict([ ('namall','Cline'),
            #        ('fixed','F'),
            # ]),

            ## Main 6.2 (1)
            dict([ ('name','Cband'+str(labB.index('Main 6.2 (1)')+1)),
                   ('fixed','F'), ]),
            dict([ ('name','WSband'+str(labB.index('Main 6.2 (1)')+1)),
                   ('fixed','F'), ]),
            dict([ ('name','WLband'+str(labB.index('Main 6.2 (1)')+1)),
                   ('fixed','F'), ]),
            ## Main 7.7 (1)
            dict([ ('name','Cband'+str(labB.index('Main 7.7 (1)')+1)),
                   ('fixed','F'), ]),
            dict([ ('name','WSband'+str(labB.index('Main 7.7 (1)')+1)),
                   ('fixed','F'), ]),
            dict([ ('name','WLband'+str(labB.index('Main 7.7 (1)')+1)),
                   ('fixed','F'), ]),
            ## Main 8.6
            dict([ ('name','Cband'+str(labB.index('Main 8.6')+1)),
                   ('fixed','F'), ]),
            dict([ ('name','WSband'+str(labB.index('Main 8.6')+1)),
                   ('fixed','F'), ]),
            dict([ ('name','WLband'+str(labB.index('Main 8.6')+1)),
                   ('fixed','F'), ]),
            ## Main 11.2
            dict([ ('name','Cband'+str(labB.index('Main 11.2')+1)),
                   ('fixed','F'), ]),
            dict([ ('name','WSband'+str(labB.index('Main 11.2')+1)),
                   ('fixed','F'), ]),
            dict([ ('name','WLband'+str(labB.index('Main 11.2')+1)),
                   ('fixed','F'), ]),

            # dict([ ('namall','lnAv'),
            #        ('value',0.),
            # ]), # LOG( 1 mag )

            dict() ]

## parinfo
##---------
Ncont = len(labQ)
Nline = len(labL)
Nband = len(labB)
Npabs = 1
Nstar = 1
Npar = 2*Ncont + 3*Nline + 4*Nband + 2

## Set indL (index for labL)
indL = []
for lab in labL:
    for tabL in TABLine:
        if tabL['label']==lab.rstrip():
            indL.append(TABLine.index(tabL))
## Set indB (index for labB)
indB = []
for lab in labB:
    for tabB in TABand:
        if tabB['label']==lab.rstrip():
            indB.append(TABand.index(tabB))

## Default setting
name = np.empty((Npar,), dtype=('<U30')) # str length <30
namall = np.empty((Npar,), dtype=('<U30')) # specific for partuning
comp = np.empty((Npar,), dtype=('<U30'))
fixed = np.array(['T' for i in range(Npar)]) # fixed by default
limited = np.array([('F','F') for i in range(Npar)])
limits = np.array([(0.,0.) for i in range(Npar)])
model = np.array(['T' for i in range(Npar)])
hyper = np.array(['F' for i in range(Npar)])
tied = np.empty((Npar,), dtype=('<U30'))
value = np.array([0. for i in range(Npar)])

## Param assignment
i0 = 0
for i in range(Ncont):
    name[i0+2*i] = 'lnMovd2'+str(i+1)
    namall[i0+2*i] = 'lnMovd2'
    value[i0+2*i] = -4. # 1.83e-2 [Msun/pc2]
    name[i0+2*i+1] = 'lnT'+str(i+1)
    namall[i0+2*i+1] = 'lnT'
    value[i0+2*i+1] = 4. # 54.60 [K]
for i in range(2*Ncont):
    comp[i0+i] = 'CONT'
i0 += 2*Ncont
for i in range(Nline):
    name[i0+3*i] = 'lnIline'+str(i+1)
    namall[i0+3*i] = 'lnIline'
    value[i0+3*i] = -20. # 2.06e-9 [MJyovsr]
    name[i0+3*i+1] = 'Cline'+str(i+1)
    namall[i0+3*i+1] = 'Cline'
    value[i0+3*i+1] = TABLine[indL[i]]['wave']
    name[i0+3*i+2] = 'Wline'+str(i+1)
    namall[i0+3*i+2] = 'Wline'
    for r in res:
        if value[i0+3*i+1]<5.:
            if r['name']=='AKARI_Ns':
                value[i0+3*i+2] = r['dwovw'] * value[i0+3*i+1]
        elif value[i0+3*i+1]<12.:
            if r['name']=='SL':
                value[i0+3*i+2] = r['dwovw'] * value[i0+3*i+1]
        else:
            if r['name']=='LL':
                value[i0+3*i+2] = r['dwovw'] * value[i0+3*i+1]
for i in range(3*Nline):
    comp[i0+i] = 'LINE'
i0 += 3*Nline
for i in range(Nband):
    name[i0+4*i] = 'lnIband'+str(i+1)
    namall[i0+4*i] = 'lnIband'
    value[i0+4*i] = -20. # 2.06e-9 [MJyovsr]
    name[i0+4*i+1] = 'Cband'+str(i+1)
    namall[i0+4*i+1] = 'Cband'
    value[i0+4*i+1] = TABand[indB[i]]['wave']
    name[i0+4*i+2] = 'WSband'+str(i+1)
    namall[i0+4*i+2] = 'WSband'
    value[i0+4*i+2] = TABand[indB[i]]['sigmaS']
    name[i0+4*i+3] = 'WLband'+str(i+1)
    namall[i0+4*i+3] = 'WLband'
    value[i0+4*i+3] = TABand[indB[i]]['sigmaL']
for i in range(4*Nband):
    comp[i0+i] = 'BAND'
i0 += 4*Nband
for i in range(Npabs):
    name[i0+i] = 'lnAv'+str(i+1)
    namall[i0+i] = 'lnAv'
    value[i0+i] = 0. # 1 [mag] ?
    comp[i0+i] = 'PABS'
i0 += Npabs
for i in range(Nstar):
    name[i0+i] = 'lnFstar'+str(i+1)
    namall[i0+i] = 'lnFstar'
    value[i0+i] = -7. # 9.12e-4 [Lsun/pc2]
    comp[i0+i] = 'STAR'

## Param tuning
partuning(dictune, Ncont, Nline, Nband,
          name, fixed, limited, limits, model, hyper, tied, value)

## Write HDF5
##------------
write_hdf5(h5_model, 'label cont', labQ, verbose=noisy)
write_hdf5(h5_model, 'label band', labB, append=True, verbose=noisy)
write_hdf5(h5_model, 'label line', labL, append=True, verbose=noisy)
write_hdf5(h5_model, 'parinfo name', name, append=True, verbose=noisy)
write_hdf5(h5_model, 'parinfo comp', comp, append=True, verbose=noisy)
write_hdf5(h5_model, 'parinfo fixed', fixed, append=True, verbose=noisy)
write_hdf5(h5_model, 'parinfo limited', limited, append=True, verbose=noisy)
write_hdf5(h5_model, 'parinfo limits', limits, append=True, verbose=noisy)
write_hdf5(h5_model, 'parinfo model', model, append=True, verbose=noisy)
write_hdf5(h5_model, 'parinfo hyper', hyper, append=True, verbose=noisy)
write_hdf5(h5_model, 'parinfo tied', tied, append=True, verbose=noisy)
write_hdf5(h5_model, 'parinfo value', value, append=True, verbose=noisy)

##--------------------------------
## Write input_fitMIR_extra.h5
##--------------------------------
Nextra = 0

## parextinfo
##------------

## Write HDF5
##------------
write_hdf5(h5_extra, 'Nextra', [Nextra], verbose=noisy)

##-----------------------------
## Write observations_fitMIR.h5 (part 2)
##-----------------------------

## Write HDF5
##------------

## These init param are supposed to be the default param in the fitting model,
## with the possibility of reasonable modifications by this script.
write_hdf5(h5_obs, 'Initial parameter label', name, append=True, verbose=noisy)
val3 = np.repeat(value[:,np.newaxis], data.shape[2], axis=1) # expand Nx
val3 = np.repeat(val3[:,:,np.newaxis], data.shape[1], axis=2) # expand Ny
write_hdf5(h5_obs, 'Initial parameter value', val3, append=True, verbose=noisy)

##--------------------------------
## Write input_fitMIR_analysis.h5
##--------------------------------

## Write HDF5
##------------

max_stars_repo_name: kxxdhdn/MISSILE

max_issues_repo_name: kxxdhdn/MISSILE

max_forks_repo_name: kxxdhdn/MISSILE


text:
"""
This script uses the cplot function in mpmath to plot the Mandelbrot set.
By default, the fp context is used for speed. The mp context could be used
to improve accuracy at extremely high zoom levels.
"""

import mpmath
import cmath

ctx = mpmath.fp
# ctx = mpmath.mp

ITERATIONS = 50
POINTS = 100000
ESCAPE_RADIUS = 8

# Full plot
RE = [-2.5, 1.5]
IM = [-1.5, 1.5]

# A pretty subplot
#RE = [-0.96, -0.80]
#IM = [-0.35, -0.2]

def mandelbrot(z):
    c = z
    for i in xrange(ITERATIONS):
        zprev = z
        z = z*z + c
        if abs(z) > ESCAPE_RADIUS:
            return ctx.exp(1j*(i + 1 - ctx.log(ctx.log(abs(z)))/ctx.log(2)))
    return 0

try:
    import psyco
    psyco.full()
except ImportError:
    pass

ctx.cplot(mandelbrot, RE, IM, points=POINTS, verbose=1)

max_stars_repo_name: upiterbarg/mpmath

max_issues_repo_name: keszybz/mpmath

max_forks_repo_name: keszybz/mpmath

reason: import mpmath


text:
"""
CKConv Implementation: taken directly from https://github.com/dwromero/ckconv
"""

from typing import Tuple, Optional

import numpy as np
import torch
import torch.fft
import torch.fft
import torch.nn as nn
import torch.nn.functional as f
from torch.nn.utils import weight_norm

def Linear1d(
        in_channels: int,
        out_channels: int,
        stride: int = 1,
        bias: bool = True,
) -> torch.nn.Module:
    """
    Implements a Linear Layer in terms of a point-wise convolution.
    """
    return nn.Conv1d(in_channels, out_channels, kernel_size=1, stride=stride, bias=bias)

def Linear2d(
        in_channels: int,
        out_channels: int,
        stride: int = 1,
        bias: bool = True,
) -> torch.nn.Module:
    """
    Implements a Linear Layer in terms of a point-wise convolution.
    """
    return nn.Conv2d(in_channels, out_channels, kernel_size=1, stride=stride, bias=bias)

def Swish():
    """
    out = x * sigmoid(x)
    """
    return Expression(lambda x: x * torch.sigmoid(x))

def Sine():
    """
    out = sin(x)
    """
    return Expression(lambda x: torch.sin(x))

class LayerNorm(nn.Module):
    def __init__(
            self,
            num_channels: int,
            eps: float = 1e-12,
    ):
        """Uses GroupNorm implementation with group=1 for speed."""
        super().__init__()
        # we use GroupNorm to implement this efficiently and fast.
        self.layer_norm = torch.nn.GroupNorm(1, num_channels=num_channels, eps=eps)

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

# From LieConv
class Expression(torch.nn.Module):
    def __init__(self, func):
        """
        Creates a torch.nn.Module that applies the function func.
        :param func: lambda function
        """
        super().__init__()
        self.func = func

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

def Multiply(
        omega_0: float,
):
    """
    out = omega_0 * x
    """
    return Expression(lambda x: omega_0 * x)

def causal_padding(
        x: torch.Tensor,
        kernel: torch.Tensor,
) -> Tuple[torch.Tensor, torch.Tensor]:
    # 1. Pad the input signal & kernel tensors.
    # Check if sizes are odd. If not, add a pad of zero to make them odd.
    if kernel.shape[-1] % 2 == 0:
        kernel = f.pad(kernel, [1, 0], value=0.0)
        # x = torch.nn.functional.pad(x, [1, 0], value=0.0)
    # 2. Perform padding on the input so that output equals input in length
    x = f.pad(x, [kernel.shape[-1] - 1, 0], value=0.0)

    return x, kernel

def causal_conv(
        x: torch.Tensor,
        kernel: torch.Tensor,
        bias: Optional[torch.Tensor] = None,
) -> torch.Tensor:
    """
    Args:
        x: (Tensor) Input tensor to be convolved with the kernel.
        kernel: (Tensor) Convolution kernel.
        bias: (Optional, Tensor) Bias tensor to add to the output.
        padding: (int) Number of zero samples to pad the input on the last dimension.
    Returns:
        (Tensor) Convolved tensor
    """

    x, kernel = causal_padding(x, kernel)
    return torch.nn.functional.conv1d(x, kernel, bias=bias, padding=0)

def causal_fftconv(
        x: torch.Tensor,
        kernel: torch.Tensor,
        bias: Optional[torch.Tensor] = None,
        double_precision: bool = False,
) -> torch.Tensor:
    """
    Args:
        x: (Tensor) Input tensor to be convolved with the kernel.
        kernel: (Tensor) Convolution kernel.
        bias: (Optional, Tensor) Bias tensor to add to the output.
        padding: (int) Number of zero samples to pad the input on the last dimension.
    Returns:
        (Tensor) Convolved tensor
    """

    x_shape = x.shape
    # 1. Handle padding of the input and the kernel to make them odd.
    x, kernel = causal_padding(x, kernel)

    # 2. Pad the kernel tensor to make them equally big. Required for fft.
    kernel = f.pad(kernel, [0, x.size(-1) - kernel.size(-1)])

    # 3. Perform fourier transform
    if double_precision:
        # We can make usage of double precision to make more accurate approximations of the convolution response.
        x = x.double()
        kernel = kernel.double()

    x_fr = torch.fft.rfft(x, dim=-1)
    kernel_fr = torch.fft.rfft(kernel, dim=-1)

    # 4. Multiply the transformed matrices:
    # (Input * Conj(Kernel)) = Correlation(Input, Kernel)
    kernel_fr = torch.conj(kernel_fr)
    output_fr = (x_fr.unsqueeze(1) * kernel_fr.unsqueeze(0)).sum(
        2
    )  # 'ab..., cb... -> ac...'

    # 5. Compute inverse FFT, and remove extra padded values
    # Once we are back in the spatial domain, we can go back to float precision, if double used.
    out = torch.fft.irfft(output_fr, dim=-1).float()

    out = out[:, :, : x_shape[-1]]

    # 6. Optionally, add a bias term before returning.
    if bias is not None:
        out = out + bias.view(1, -1, 1)

    return out

class KernelNet(torch.nn.Module):
    def __init__(
            self,
            in_channels: int,
            out_channels: int,
            hidden_channels: int,
            activation_function: str,
            norm_type: str,
            dim_linear: int,
            bias: bool,
            omega_0: float,
            weight_dropout: float,
    ):
        """
        Creates an 3-layer MLP, which parameterizes a convolutional kernel as:
        relative positions -> hidden_channels -> hidden_channels -> in_channels * out_channels
        :param in_channels:  Dimensionality of the relative positions (Default: 1).
        :param out_channels:  input channels * output channels of the resulting convolutional kernel.
        :param hidden_channels: Number of hidden units.
        :param activation_function: Activation function used.
        :param norm_type: Normalization type used.
        :param dim_linear:  Spatial dimension of the input, e.g., for audio = 1, images = 2 (only 1 suported).
        :param bias:  If True, adds a learnable bias to the layers.
        :param omega_0: Value of the omega_0 value (only used in Sine networks).
        :param weight_dropout: Dropout rate applied to the sampled convolutional kernel.
        """
        super().__init__()

        is_siren = activation_function == "Sine"
        w_dp = weight_dropout != 0.0

        Norm = {
            "BatchNorm": torch.nn.BatchNorm1d,
            "LayerNorm": LayerNorm,
            "": torch.nn.Identity,
        }[norm_type]
        ActivationFunction = {
            "ReLU": torch.nn.ReLU,
            "LeakyReLU": torch.nn.LeakyReLU,
            "Swish": Swish,
            "Sine": Sine,
        }[activation_function]
        Linear = {1: Linear1d, 2: Linear2d}[dim_linear]

        self.kernel_net = torch.nn.Sequential(
            weight_norm(Linear(in_channels, hidden_channels, bias=bias)),
            Multiply(omega_0) if is_siren else torch.nn.Identity(),
            Norm(hidden_channels) if not is_siren else torch.nn.Identity(),
            ActivationFunction(),
            weight_norm(Linear(hidden_channels, hidden_channels, bias=bias)),
            Multiply(omega_0) if is_siren else torch.nn.Identity(),
            Norm(hidden_channels) if not is_siren else torch.nn.Identity(),
            ActivationFunction(),
            weight_norm(Linear(hidden_channels, out_channels, bias=bias)),
            torch.nn.Dropout(p=weight_dropout) if w_dp else torch.nn.Identity(),
        )

        # initialize the kernel function
        self.initialize(
            mean=0.0,
            variance=0.01,
            bias_value=0.0,
            is_siren=(activation_function == "Sine"),
            omega_0=omega_0,
        )

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

    def initialize(self, mean, variance, bias_value, is_siren, omega_0):

        if is_siren:
            # Initialization of SIRENs
            net_layer = 1
            for (i, m) in enumerate(self.modules()):
                if (
                        isinstance(m, torch.nn.Conv1d)
                        or isinstance(m, torch.nn.Conv2d)
                        or isinstance(m, torch.nn.Linear)
                ):
                    if net_layer == 1:
                        m.weight.data.uniform_(
                            -1, 1
                        )  # Normally (-1, 1) / in_dim but we only use 1D inputs.
                        # Important! Bias is not defined in original SIREN implementation!
                        net_layer += 1
                    else:
                        m.weight.data.uniform_(
                            -np.sqrt(6.0 / m.weight.shape[1]) / omega_0,
                            # the in_size is dim 2 in the weights of Linear and Conv layers
                            np.sqrt(6.0 / m.weight.shape[1]) / omega_0,
                        )
                    # Important! Bias is not defined in original SIREN implementation
                    if m.bias is not None:
                        m.bias.data.uniform_(-1.0, 1.0)
        else:
            # Initialization of ReLUs
            net_layer = 1
            intermediate_response = None
            for (i, m) in enumerate(self.modules()):
                if (
                        isinstance(m, torch.nn.Conv1d)
                        or isinstance(m, torch.nn.Conv2d)
                        or isinstance(m, torch.nn.Linear)
                ):
                    m.weight.data.normal_(
                        mean,
                        variance,
                    )
                    if m.bias is not None:

                        if net_layer == 1:
                            # m.bias.data.fill_(bias_value)
                            range = torch.linspace(-1.0, 1.0, steps=m.weight.shape[0])
                            bias = -range * m.weight.data.clone().squeeze()
                            m.bias = torch.nn.Parameter(bias)

                            intermediate_response = [
                                m.weight.data.clone(),
                                m.bias.data.clone(),
                            ]
                            net_layer += 1

                        elif net_layer == 2:
                            range = torch.linspace(-1.0, 1.0, steps=m.weight.shape[0])
                            range = range + (range[1] - range[0])
                            range = (
                                    range * intermediate_response[0].squeeze()
                                    + intermediate_response[1]
                            )

                            bias = -torch.einsum(
                                "oi, i -> o", m.weight.data.clone().squeeze(), range
                            )
                            m.bias = torch.nn.Parameter(bias)

                            net_layer += 1

                        else:
                            m.bias.data.fill_(bias_value)

class CKConv(torch.nn.Module):
    def __init__(
            self,
            in_channels: int,
            out_channels: int,
            hidden_channels: int,
            activation_function: str,
            norm_type: str,
            dim_linear: int,
            bias: bool,
            omega_0: float,
            weight_dropout: float,
    ):
        """
        Creates a Continuous Kernel Convolution.
        :param in_channels: Number of channels in the input signal
        :param out_channels: Number of channels produced by the convolution
        :param hidden_channels: Number of hidden units in the network parameterizing the ConvKernel (KernelNet).
        :param activation_function: Activation function used in KernelNet.
        :param norm_type: Normalization type used in KernelNet. (only for non-Sine KernelNets).
        :param dim_linear: patial dimension of the input, e.g., for audio = 1, images = 2 (only 1 suported).
        :param bias: If True, adds a learnable bias to the output.
        :param omega_0: Value of the omega_0 value of the KernelNet. (only for non-Sine KernelNets).
        :param weight_dropout: Dropout rate applied to the sampled convolutional kernels.
        """
        super().__init__()
        self.Kernel = KernelNet(
            dim_linear,
            out_channels * in_channels,
            hidden_channels,
            activation_function,
            norm_type,
            dim_linear,
            bias,
            omega_0,
            weight_dropout,
        )

        if bias:
            self.bias = torch.nn.Parameter(torch.Tensor(out_channels))
            self.bias.data.fill_(value=0.0)
        else:
            self.bias = None

        # Non-persistent values
        self.rel_positions = None
        self.sigma = None
        self.sr_change = 1.0

        self.register_buffer("train_length", torch.zeros(1).int(), persistent=True)
        self.register_buffer("conv_kernel", torch.zeros(in_channels), persistent=False)

    def forward(self, x):
        # Construct kernel
        x_shape = x.shape

        rel_pos = self.handle_rel_positions(x)
        conv_kernel = self.Kernel(rel_pos).view(-1, x_shape[1], *x_shape[2:])

        # ---- Different samling rate --------
        # If freq test > freq test, smooth out high-freq elements.
        if self.sigma is not None:
            from math import pi, sqrt, exp

            n = int(1 / self.sr_change) * 2 + 1
            h = n // 2
            G = (
                lambda x: 1
                          / (self.sigma * sqrt(2 * pi))
                          * exp(-float(x) ** 2 / (2 * self.sigma ** 2))
            )

            smoothing_ker = [G(x) for x in range(-h, h + 1)]
            smoothing_ker = torch.Tensor(smoothing_ker).cuda().unsqueeze(0).unsqueeze(0)
            conv_kernel[:, :, h:-h] = torch.conv1d(
                conv_kernel.view(-1, 1, *x_shape[2:]), smoothing_ker, padding=0
            ).view(*conv_kernel.shape[:-1], -1)
        # multiply by the sr_train / sr_test
        if self.sr_change != 1.0:
            conv_kernel *= self.sr_change
        # ------------------------------------

        # For computation of "weight_decay"
        self.conv_kernel = conv_kernel

        # We have noticed that the results of fftconv become very noisy when the length of
        # the input is very small ( < 50 samples). As this might occur when we use subsampling,
        # we replace causal_fftconv by causal_conv in settings where this occurs.
        if x_shape[-1] < self.train_length.item():
            # Use spatial convolution:
            return causal_conv(x, conv_kernel, self.bias)
        else:
            # Otherwise use fft convolution:
            return causal_fftconv(x, conv_kernel, self.bias)

    def handle_rel_positions(self, x):
        """
        Handles the vector or relative positions which is given to KernelNet.
        """
        if self.rel_positions is None:
            if self.train_length[0] == 0:
                # The ckconv has not been trained yet. Set maximum length to be 1.
                self.train_length[0] = x.shape[-1]

            # Calculate the maximum relative position based on the length of the train set,
            # and the current length of the input.
            max_relative_pos = self.calculate_max(
                self.train_length.item(), current_length=x.shape[-1]
            )

            # Creates the vector of relative positions.
            self.rel_positions = (
                torch.linspace(-1.0, max_relative_pos, x.shape[-1])
                    .cuda()
                    .unsqueeze(0)
                    .unsqueeze(0)
            )  # -> With form: [batch_size=1, in_channels=1, x_dimension]

            # calculate and save the sr ratio for later
            if self.train_length.item() > x.shape[-1]:
                self.sr_change = round(self.train_length.item() / x.shape[-1])
            else:
                self.sr_change = 1 / round(x.shape[-1] / self.train_length.item())

            # if new signal has higher frequency
            if self.sr_change < 1:
                self.sigma = 0.5

        return self.rel_positions

    @staticmethod
    def calculate_max(
            train_length: int,
            current_length: int,
    ) -> float:
        """
        Calculates the maximum relative position for the current length based on the input length.
        This is used to avoid kernel misalignment (see Appx. D.2).
        :param train_length: Input length during training.
        :param current_length: Current input length.
        :return: Returns the max relative position for the calculation of the relative
                 positions vector. The max. of train is always equal to 1.
        """
        # get sampling rate ratio
        if train_length > current_length:
            sr_change = round(train_length / current_length)
        else:
            sr_change = 1 / round(current_length / train_length)

        # get step sizes (The third parameter of torch.linspace).
        train_step = 2.0 / (train_length - 1)
        current_step = train_step * sr_change

        # Calculate the maximum relative position.
        if sr_change > 1:
            substract = (train_length - 1) % sr_change
            max_relative_pos = 1 - substract * train_step
        else:
            add = (current_length - 1) % (1 / sr_change)
            max_relative_pos = 1 + add * current_step
        return max_relative_pos

class CKBlock(torch.nn.Module):
    def __init__(
            self,
            in_channels: int,
            out_channels: int,
            kernelnet_hidden_channels: int,
            kernelnet_activation_function: str,
            kernelnet_norm_type: str,
            dim_linear: int,
            bias: bool,
            omega_0: bool,
            dropout: float,
            weight_dropout: float,
    ):
        """
        Creates a Residual Block with CKConvs as:
        ( Follows the Residual Block of Bai et. al., 2017 )
        input
         | ---------------|
         CKConv           |
         LayerNorm        |
         ReLU             |
         DropOut          |
         |                |
         CKConv           |
         LayerNorm        |
         ReLU             |
         DropOut          |
         + <--------------|
         |
         ReLU
         |
         output
        :param in_channels:  Number of channels in the input signal
        :param out_channels:  Number of output (and hidden) channels of the block.
        :param kernelnet_hidden_channels: Number of hidden units in the KernelNets of the CKConvs.
        :param kernelnet_activation_function: Activation function used in the KernelNets of the CKConvs.
        :param kernelnet_norm_type: Normalization type used in the KernelNets of the CKConvs (only for non-Sine KernelNets).
        :param dim_linear:  Spatial dimension of the input, e.g., for audio = 1, images = 2 (only 1 suported).
        :param bias:  If True, adds a learnable bias to the output.
        :param omega_0: Value of the omega_0 value of the KernelNets. (only for non-Sine KernelNets).
        :param dropout: Dropout rate of the block
        :param weight_dropout: Dropout rate applied to the sampled convolutional kernels.
        """
        super().__init__()

        # CKConv layers
        self.cconv1 = CKConv(
            in_channels,
            out_channels,
            kernelnet_hidden_channels,
            kernelnet_activation_function,
            kernelnet_norm_type,
            dim_linear,
            bias,
            omega_0,
            weight_dropout,
        )
        self.cconv2 = CKConv(
            out_channels,
            out_channels,
            kernelnet_hidden_channels,
            kernelnet_activation_function,
            kernelnet_norm_type,
            dim_linear,
            bias,
            omega_0,
            weight_dropout,
        )
        # Norm layers
        self.norm1 = LayerNorm(out_channels)
        self.norm2 = LayerNorm(out_channels)

        # Dropout
        self.dp = torch.nn.Dropout(dropout)

        shortcut = []
        if in_channels != out_channels:
            shortcut.append(Linear1d(in_channels, out_channels))
        self.shortcut = torch.nn.Sequential(*shortcut)

    def forward(self, x):
        shortcut = self.shortcut(x)
        out = self.dp(torch.relu(self.norm1(self.cconv1(x))))
        out = torch.relu(self.dp(torch.relu(self.norm2(self.cconv2(out)))) + shortcut)
        return out

class CKCNN(torch.nn.Module):
    def __init__(
            self,
            in_channels: int,
            hidden_channels: int,
            num_blocks: int,  # 2
            kernelnet_hidden_channels: int,
            kernelnet_activation_function: str,
            kernelnet_norm_type: str,
            dim_linear: int,
            bias: bool,
            omega_0: bool,  # sensitive to this param: good values <= 70
            dropout: float,
            weight_dropout: float,
            pool: bool,  # Always False in our experiments.
    ):
        super(CKCNN, self).__init__()

        blocks = []
        for i in range(num_blocks):
            block_in_channels = in_channels if i == 0 else hidden_channels
            blocks.append(
                CKBlock(
                    block_in_channels,
                    hidden_channels,
                    kernelnet_hidden_channels,
                    kernelnet_activation_function,
                    kernelnet_norm_type,
                    dim_linear,
                    bias,
                    omega_0,
                    dropout,
                    weight_dropout,
                )
            )
            if pool:
                blocks.append(torch.nn.MaxPool1d(kernel_size=2))
        self.backbone = torch.nn.Sequential(*blocks)

    def forward(self, x, *args, **kwargs):
        # Change from (B, S, C) -> (B, C, S)
        x = x.transpose(1, 2)
        return self.backbone(x)

class CopyMemory_CKCNN(CKCNN):
    def __init__(
            self,
            in_channels: int,
            hidden_channels: int,
            num_blocks: int,
            kernelnet_hidden_channels: int,
            kernelnet_activation_function: str,
            kernelnet_norm_type: str,
            dim_linear: int,
            bias: bool,
            omega_0: bool,
            dropout: float,
            weight_dropout: float,
            pool: bool,
    ):
        super().__init__(
            in_channels,
            hidden_channels,
            num_blocks,
            kernelnet_hidden_channels,
            kernelnet_activation_function,
            kernelnet_norm_type,
            dim_linear,
            bias,
            omega_0,
            dropout,
            weight_dropout,
            pool,
        )

        self.finallyr = torch.nn.Linear(in_features=hidden_channels, out_features=10)
        # Initialize finallyr
        self.finallyr.weight.data.normal_(
            mean=0.0,
            std=0.01,
        )
        self.finallyr.bias.data.fill_(value=0.0)

    def forward(self, x, *args, **kwargs):
        # Change from (B, S, C) -> (B, C, S)
        x = x.transpose(1, 2)
        out = self.backbone(x)
        out = self.finallyr(out.transpose(1, 2))
        return out

class AddProblem_CKCNN(CKCNN):
    def __init__(
            self,
            in_channels: int,
            hidden_channels: int,
            num_blocks: int,
            kernelnet_hidden_channels: int,
            kernelnet_activation_function: str,
            kernelnet_norm_type: str,
            dim_linear: int,
            bias: bool,
            omega_0: bool,
            dropout: float,
            weight_dropout: float,
            pool: bool,
    ):
        super().__init__(
            in_channels,
            hidden_channels,
            num_blocks,
            kernelnet_hidden_channels,
            kernelnet_activation_function,
            kernelnet_norm_type,
            dim_linear,
            bias,
            omega_0,
            dropout,
            weight_dropout,
            pool,
        )

        self.finallyr = torch.nn.Linear(in_features=hidden_channels, out_features=1)
        # Initialize finallyr
        self.finallyr.weight.data.normal_(
            mean=0.0,
            std=0.01,
        )
        self.finallyr.bias.data.fill_(value=0.0)

    def forward(self, x, *args, **kwargs):
        # Change from (B, S, C) -> (B, C, S)
        x = x.transpose(1, 2)
        out = self.backbone(x)
        out = self.finallyr(out[:, :, -1])
        return out

class ClassificationCKCNN(CKCNN):
    def __init__(
            self,
            d_input: int,
            d_output: int,
            hidden_channels: int,
            num_blocks: int,
            kernelnet_hidden_channels: int,
            kernelnet_activation_function: str,
            kernelnet_norm_type: str,
            dim_linear: int,
            bias: bool,
            omega_0: bool,
            dropout: float,
            weight_dropout: float,
            pool: bool,
            wd: float,
            **kwargs,
    ):
        super().__init__(
            d_input,
            hidden_channels,
            num_blocks,
            kernelnet_hidden_channels,
            kernelnet_activation_function,
            kernelnet_norm_type,
            dim_linear,
            bias,
            omega_0,
            dropout,
            weight_dropout,
            pool,
        )
        if 'l_output' in kwargs and kwargs['l_output'] > 1:
            d_output = kwargs['l_output']

        self.finallyr = torch.nn.Linear(
            in_features=hidden_channels, out_features=d_output
        )
        # Initialize finallyr
        self.finallyr.weight.data.normal_(
            mean=0.0,
            std=0.01,
        )
        self.finallyr.bias.data.fill_(value=0.0)

        self.wd = LnLoss(wd, 2)

    def forward(self, x, *args, **kwargs):
        # Change from (B, S, C) -> (B, C, S)
        x = x.transpose(1, 2)
        out = self.backbone(x)
        out = self.finallyr(out[:, :, -1])
        return out

    # def decay(self):
    def loss(self):
        return self.wd.forward(model=self)

class LnLoss(torch.nn.Module):
    def __init__(
            self,
            weight_loss: float,
            norm_type: int,
    ):
        """
        Computes the Ln loss on the CKConv kernels in a CKCNN.
        :param weight_loss: Specifies the weight with which the loss will be summed to the total loss.
        :param norm_type: Type of norm, e.g., 1 = L1 loss, 2 = L2 loss, ...
        """
        super(LnLoss, self).__init__()
        self.weight_loss = weight_loss
        self.norm_type = norm_type

    def forward(
            self,
            model: CKConv,
    ):
        loss = 0.0
        # Go through modules that are instances of CKConvs and gather the sampled filters
        for m in model.modules():
            if not isinstance(m, CKConv):
                continue
            loss = loss + m.conv_kernel.norm(self.norm_type)
            loss = loss + m.bias.norm(self.norm_type)

        loss = self.weight_loss * loss
        return loss

max_stars_repo_name: dumpmemory/state-spaces

max_issues_repo_name: dumpmemory/state-spaces

max_forks_repo_name: MikeOwino/state-spaces


text:
'''
USAGE:
python preprocess_image.py --num-images 1200
'''
import os
import cv2
import random
import numpy as np
import argparse
from tqdm import tqdm
parser = argparse.ArgumentParser()
parser.add_argument('-n', '--num-images', default=1200, type=int,
    help='number of images to preprocess for each category')
args = vars(parser.parse_args())
print(f"Preprocessing {args['num_images']} from each category...")

# get all the directory paths
dir_paths = os.listdir('asl_alphabet_train/asl_alphabet_train')
dir_paths.sort()
root_path = 'asl_alphabet_train/asl_alphabet_train'

# get --num-images images from each category
for idx, dir_path  in tqdm(enumerate(dir_paths), total=len(dir_paths)):
    all_images = os.listdir(f"{root_path}/{dir_path}")
    os.makedirs(f"preprocessed_image/{dir_path}", exist_ok=True)
    for i in range(args['num_images']): # how many images to preprocess for each category
        # generate a random id between 0 and 2999
        rand_id = (random.randint(0, 2999))
        image = cv2.imread(f"{root_path}/{dir_path}/{all_images[rand_id]}")
        image = cv2.resize(image, (224, 224))
        cv2.imwrite(f"preprocessed_image/{dir_path}/{dir_path}{i}.jpg", image)
print('DONE')

max_stars_repo_name: MhmudAlpurd/IC-pytorchlite

max_issues_repo_name: MhmudAlpurd/IC-pytorchlite

max_forks_repo_name: MhmudAlpurd/IC-pytorchlite


text:
import cPickle
import gzip
import os

import astropy.units as u
from astropy.coordinates import SkyCoord

class ChandraSourceCatalog(object):
    def __init__(self):

        # Find the catalog data file

        # TODO: Put it in the region repository

        catalog_filename = 'chandra_csc_1.1.pickle.gz'

        # Define the search paths
        search_paths = ['.', os.path.dirname(__file__)]

        catalog_file = None

        for path in search_paths:

            catalog_file = os.path.abspath(os.path.join(path, catalog_filename))

            if os.path.exists(catalog_file):
                # Found it
                break

        assert catalog_file is not None, "Could not find catalog file in %s" % ",".join(search_paths)

        # Read the chandra source catalog from the CSV file

        f = gzip.GzipFile(catalog_file)

        data = cPickle.load(f)

        self._catalog = data['data_frame']
        self._sky_coords = data['sky_coords']

    def _compute_distances(self, ra, dec):

        # Instance the SkyCoord instance

        cone_center = SkyCoord(ra=ra, dec=dec, unit='deg')

        # Find all sources within the requested cone

        distances = self._sky_coords.separation(cone_center).to(u.arcmin)

        return distances

    def cone_search(self, ra, dec, radius, unit='arcmin'):
        """
        Find sources within the given radius of the given position

        :param ra: R.A. of position
        :param dec: Dec of position
        :param radius: radius
        :param unit: units to use for the radius (default: arcmin)
        :return: a pandas DataFrame containing the sources within the given radius
        """

        distances = self._compute_distances(ra, dec)  # arcmin

        idx = distances <= radius * u.Unit(unit)

        # Return the results

        # Copy the array, to avoid returning a slice instead

        results = self._catalog.copy().loc[idx]

        # Add the distance column
        results['distance'] = distances[idx]

        return results

    def find_closest_source(self, ra, dec):
        """
        Finds the closest source to the given position

        :param ra:
        :param dec:
        :return:
        """

        distances = self._compute_distances(ra, dec)  # arcmin

        temp_catalog = self._catalog.copy()

        temp_catalog['distance'] = distances

        src_id = temp_catalog['distance'].argmin()

        return temp_catalog.loc[src_id, :]

    def find_variable_sources(self, ra, dec, radius, unit='arcmin', column='var_flag'):
        """
        Find all variable sources within the given cone

        :param ra: R.A. of position
        :param dec: Dec of position
        :param radius: radius
        :param unit: units to use for the radius (default: arcmin)
        :return: a pandas DataFrame containing the variable sources within the given radius
        """

        # Get all sources within the cone

        temp_results = self.cone_search(ra, dec, radius, unit=unit)

        # Now select only the variable sources
        idx = temp_results[column] == True

        # Get a copy
        results = temp_results.copy().loc[idx]

        return results

    def find_closest_variable_source(self, ra, dec, column='var_flag'):
        """
        Finds the closest source to the given position

        :param ra:
        :param dec:
        :return:
        """

        distances = self._compute_distances(ra, dec)  # arcmin

        temp_catalog = self._catalog.copy()

        temp_catalog['distance'] = distances

        # Select only variable sources
        idx = temp_catalog[column] == True

        variable_sources = temp_catalog.copy().loc[idx]

        src_id = variable_sources['distance'].argmin()

        return variable_sources.loc[src_id, :]

reason: import astropy,from astropy


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

'''
@Author:        ZM
@Date and Time: 2020/12/16 6:52
@File:          generator.py
'''

import numpy as np

def generator(dataset, batch_size=64, shuffle=True, drop_last=False):
    text, label = dataset[0]
    text_size, num_classes = text.shape, label.shape

    true_examples = len(dataset)
    rd_index = np.arange(true_examples)
    false_examples = true_examples // batch_size * batch_size
    remain_examples = true_examples - false_examples

    i = 0
    while True:
        real_batch_size = batch_size
        if remain_examples != 0 and drop_last is False and i == false_examples:
            real_batch_size = remain_examples

        batch_texts = np.empty((real_batch_size, *text_size), dtype='float32')
        batch_labels = np.empty((real_batch_size, *num_classes), dtype='float32')

        for b in range(real_batch_size):
            if shuffle and i == 0:
                np.random.shuffle(rd_index)
                dataset.texts = dataset.texts[rd_index]
                dataset.labels = dataset.labels[rd_index]

            batch_texts[b], batch_labels[b] = dataset[i]

            if remain_examples != 0 and drop_last is True:
                i = (i + 1) % false_examples
            else:
                i = (i + 1) % true_examples

        yield batch_texts, batch_labels


text:
import matplotlib
matplotlib.use('Agg')
import pylab
import sys
import numpy
from utils import stats_utils
import matplotlib.colors as colors
import matplotlib.cm as cmx
from math import log10,ceil
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from numpy.random import randint
import pandas
import os
from parsers import parse_midas_data
import pickle

species_names = parse_midas_data.parse_good_species_list()
# load species info
files={}

for species_name in  species_names:
    if (os.path.exists(os.path.expanduser('~/tmp_intermediate_files/ortholog_%s.dat' % species_name))):
        files[species_name]=pickle.load(open(os.path.expanduser('~/tmp_intermediate_files/ortholog_%s.dat' % species_name),'rb'))
    else:
        print species_name

# iterate through the species, add the counts to vectors
num_genes_with_match_within_host=[]
num_genes_with_match_random=[]

for species_name in files:
    num_genes_with_match_within_host.append(files[species_name]['num_genes_with_match_within_host'])
    num_genes_with_match_random.append(files[species_name]['num_genes_with_match_random'])

#plot

pylab.figure()
pylab.xlabel('number of within-host genes with a 100% match')
pylab.ylabel('number of random genes with a 100% match')
pylab.title('Number of genes with a 100% match to an off-target species')

pylab.plot(num_genes_with_match_within_host, num_genes_with_match_random, "ro")
pylab.plot([0,800],[0,800] , 'k-')
pylab.savefig('%s/within_vs_random_gene_orthologs.png' % (parse_midas_data.analysis_directory), bbox_inches='tight')

reason: import numpy,from numpy


text:
#%%
import os

try:
    os.chdir('/Volumes/GoogleDrive/My Drive/python_code/connectome_tools/')
    print(os.getcwd())
except:

    pass

#%%
import sys
sys.path.append("/Volumes/GoogleDrive/My Drive/python_code/connectome_tools/")

import pandas as pd
import numpy as np
import connectome_tools.process_matrix as promat
import math
import matplotlib.pyplot as plt
import seaborn as sns
from tqdm import tqdm
import pymaid
from pymaid_creds import url, name, password, token

# convert pair-sorted brain/sensories matrix to binary matrix based on synapse threshold
matrix_ad = pd.read_csv('data/axon-dendrite.csv', header=0, index_col=0)
matrix_dd = pd.read_csv('data/dendrite-dendrite.csv', header=0, index_col=0)
matrix_aa = pd.read_csv('data/axon-axon.csv', header=0, index_col=0)
matrix_da = pd.read_csv('data/dendrite-axon.csv', header=0, index_col=0)

matrix = matrix_ad + matrix_dd + matrix_aa + matrix_da

# the columns are string by default and the indices int; now both are int
matrix_ad.columns = pd.to_numeric(matrix_ad.columns)
matrix_dd.columns = pd.to_numeric(matrix_dd.columns)
matrix_aa.columns = pd.to_numeric(matrix_aa.columns)
matrix_da.columns = pd.to_numeric(matrix_da.columns)
matrix.columns = pd.to_numeric(matrix.columns)

# import pair list CSV, manually generated
pairs = pd.read_csv('data/pairs-2020-05-08.csv', header = 0)
paired = pairs.values.flatten()

# %%
rm = pymaid.CatmaidInstance(url, token, name, password)

# pull sensory annotations and then pull associated skids
sensories = pymaid.get_annotated('mw brain inputs')

# %%

testinput = pymaid.get_skids_by_annotation("mw ORN 2nd_order PN")

# identify pairs in a-d graph

def get_paired_skids(skid, pairList):
# returns paired skids in array [left, right]; can input either left or right skid of a pair to identify

    if(skid in pairList["leftid"].values):
        pair_right = pairList["rightid"][pairList["leftid"]==skid].iloc[0]
        pair_left = skid

    if(skid in pairList["rightid"].values):
        pair_left = pairList["leftid"][pairList["rightid"]==skid].iloc[0]
        pair_right = skid

    if((skid in pairList["leftid"].values) == False and (skid in pairList["rightid"].values) == False):
        print("skid %i is not in paired list" % (skid))
        return(0)

    return([pair_left, pair_right])

def extract_pairs_from_list(skids, pairList):
    pairs = []
    for i in skids:
        if(int(i) in pairList["leftid"].values):
            pair = get_paired_skids(int(i), pairList)
            pairs.append({'leftid': pair[0], 'rightid': pair[1]})

    pairs = pd.DataFrame(pairs)
    return(pairs)

testpairs = extract_pairs_from_list(testinput, pairs)

# look for downstream pairs in a-d graph

# identify LNs and descending in downstream pairs

# repeat

# %%


text:
import pathlib
from typing import Union,Tuple
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import re
import math

class QuestionnaireAnalysis:
    """
    Reads and analyzes data generated by the questionnaire experiment.
    Should be able to accept strings and pathlib.Path objects.
    """

    def __init__(self, data_fname: Union[pathlib.Path, str]):
        self.data = None
        if type(data_fname)== str:
            self.data_fname= pathlib.Path(data_fname)
        else:
            self.data_fname= data_fname
        if not self.data_fname.is_file():
            raise  ValueError('file does not exist')

    def read_data(self):
        """Reads the json data located in self.data_fname into memory, to
        the attribute self.data.
        """
        self.data= pd.read_json(self.data_fname)

    def show_age_distrib(self):
        df= self.data
        ages = df['age'].to_numpy()
        edges = np.array([i*10 for i in range(11)])
        list_of_ages = plt.hist(ages,bins=edges)[0]
        plt.show()
        return  (list_of_ages,edges)
    def remove_rows_without_mail(self):

        df = self.data
        self.data =  df[df['email'].map(validEmail)]
        self.data = self.data.reset_index(drop=True)
        return self.data

    def fill_na_with_mean(self) -> Tuple[pd.DataFrame, np.ndarray]:
        df= self.data
        badRows = set()
        for index, row in df.iterrows():
            for i in range(1,6):
                if math.isnan(row['q{0}'.format(i)]):
                    badRows.add(index)
                    df['q{0}'.format(i)][index] = np.nanmean([row['q1'],row['q2'],row['q3'],row['q4'],row['q5']])
        return (df,np.array(list(badRows)))

    def score_subjects(self, maximal_nans_per_sub: int = 1) -> pd.DataFrame:
        df = self.data
        values = []
        for index, row in df.iterrows():
            count = 0
            for i in range(1,6):
                if math.isnan(row['q{0}'.format(i)]):
                    count +=1
            if count > maximal_nans_per_sub:
                values.append(None)
            else:
                values.append(int(np.nanmean([row['q1'],row['q2'],row['q3'],row['q4'],row['q5']])))
        df['score'] = values
        df['score'] = df['score'].astype(pd.UInt8Dtype())
        self.data = df
        return self.data

    def correlate_gender_age(self) -> pd.DataFrame:
        #index = pd.MultiIndex.from_arrays(self.data, names=('gender', 'age'))
        df = self.data[['gender', 'age','q1','q2','q3','q4','q5']]
        df.age = df.age.apply(lambda x: x > 40)
        df = df.set_index(['gender', 'age'])
        df = df.groupby(['gender', 'age']).mean()
        return df

def validEmail(email):
    regex = '.+@.+\..+'
    pattern = re.compile(regex)
    return re.match(pattern,email) is not None