text:
from typing import Callable, List, Optional, Tuple

import numpy as np

from robustnessgym.core.constants import TRANSFORMATION
from robustnessgym.core.identifier import Identifier
from robustnessgym.core.slice import SliceDataPanel as DataPanel
from robustnessgym.slicebuilders.slicebuilder import SliceBuilder

class Transformation(SliceBuilder):
    def __init__(
        self,
        num_transformed: int = None,
        identifiers: List[Identifier] = None,
        apply_fn: Callable = None,
        category: str = None,
    ):
        assert (
            num_transformed if not identifiers else True
        ), "Must pass in num_transformed if no identifiers are given."

        super(Transformation, self).__init__(
            identifiers=[
                Identifier(
                    _name=f"{self.__class__.__name__}-{i + 1}",
                )
                for i in range(num_transformed)
            ]
            if not identifiers
            else identifiers,
            category=category if category else TRANSFORMATION,
            apply_fn=apply_fn,
        )

    @property
    def num_transformed(self):
        return self.num_slices

    def apply(
        self,
        batch: DataPanel,
        columns: List[str],
        skeleton_batches: List[DataPanel],
        slice_membership: np.ndarray,
        *args,
        **kwargs,
    ) -> Tuple[List[DataPanel], np.ndarray]:
        raise NotImplementedError

    def process_batch(
        self,
        dp: DataPanel,
        columns: List[str],
        *args,
        **kwargs,
    ) -> Tuple[List[DataPanel], Optional[np.ndarray]]:

        # Determine the size of the batch
        batch_size = len(dp[list(dp.keys())[0]])

        # Construct the matrix of slice labels: (batch_size x num_slices)
        slice_membership = np.ones((batch_size, self.num_slices), dtype=np.int32)

        # Uncache the batch to construct the skeleton for transformed batches
        skeleton_batches = [
            DataPanel.uncached_batch(dp) for _ in range(self.num_slices)
        ]

        # Set the index for the skeleton batches
        for j, skeleton_batch in enumerate(skeleton_batches):
            # skeleton_batch.update(
            #     lambda x: {'index': f"{x['index']}-{self.identifiers[j]}"}
            # )
            skeleton_batch["index"] = [
                f"{idx}-{self.identifiers[j]}" for idx in skeleton_batch["index"]
            ]

        # Apply the SliceBuilder's core functionality: use positional args
        try:
            transformed_batches, slice_membership = self.apply(
                dp,
                columns,
                skeleton_batches,
                slice_membership,
                *args,
                **kwargs,
            )
        except TypeError:
            self.apply(dp, columns, *args, **kwargs)

        # Remove transformed examples where slice_membership[i, :] = 0 before returning
        transformed_batches = [
            self.filter_batch_by_slice_membership(
                batch=transformed_batch,
                slice_membership=slice_membership[:, j : j + 1],
            )[0]
            for j, transformed_batch in enumerate(transformed_batches)
        ]

        return transformed_batches, slice_membership

class SingleColumnTransformation(Transformation):
    def single_column_apply(self, column_batch: List) -> List[List]:
        return NotImplemented(
            "Implement single_column_apply to use this transformation."
        )

    def apply(
        self,
        batch: DataPanel,
        columns: List[str],
        skeleton_batches: List[DataPanel],
        slice_membership: np.ndarray,
        *args,
        **kwargs,
    ) -> Tuple[List[DataPanel], np.ndarray]:

        # Independently apply the transformation over the columns
        for column in columns:
            try:
                # Apply
                transformed_batch = self.single_column_apply(
                    column_batch=batch[column],
                )

                assert len(transformed_batch) == len(
                    batch[column]
                ), "Must output one list of augmentations per example."

                # Store the transformed text in the skeleton batches
                for i in range(slice_membership.shape[0]):
                    for j, transformed in enumerate(transformed_batch[i]):
                        skeleton_batches[j][column][i] = transformed
            except:  # noqa
                # Unable to transform: set all slice membership labels to zero
                slice_membership[:, :] = 0

        return skeleton_batches, slice_membership

max_stars_repo_name: jessevig/robustness-gym

max_issues_repo_name: jessevig/robustness-gym

max_forks_repo_name: jessevig/robustness-gym


text:
"""
Using IBIS8542Model
===================
This is an example showing how to fit an array of
spectra using the :class:`mcalf.models.IBIS8542Model`
class, and plot and export the results.
"""

#%%
# Generate sample data
# --------------------
#
# Create spectral grid
# ~~~~~~~~~~~~~~~~~~~~
#
# First, we shall generate some random sample data to
# demonstrate the API.
# The randomly generated spectra are not intended to
# be representative of observed spectra numerically,
# they just have a similar shape.

import numpy as np
import matplotlib.pyplot as plt

from mcalf.profiles.voigt import voigt, double_voigt
from mcalf.visualisation import plot_spectrum

# Create demo wavelength grid
w = np.linspace(8541, 8543, 20)

#%%
# We shall generate demo spectra in a 2x3 grid.
# Half of the spectra will have one spectral
# component (absorption only) and the other
# half will have two spectral components
# (mixed absorption and emission components).
#
# Demo spectral components will be modelled as
# Voigt functions with randomly generated
# parameters, each within a set range of values.

def v(classification, w, *args):
    """Voigt function wrapper."""
    if classification < 1.5:  # absorption only
        return voigt(w, *args[:4], args[-1])
    return double_voigt(w, *args)  # absorption + emission

def s():
    """Generate random spectra for a 2x3 grid."""
    np.random.seed(0)  # same spectra every time
    p = np.random.rand(9, 6)  # random Voigt parameters
    # 0  1  2  3  4  5  6  7  8  # p index
    # a1 b1 s1 g1 a2 b2 s2 g2 d  # Voigt parameter
    # absorption |emission   |background

    p[0] = 100 * p[0] - 1000  # absorption amplitude
    p[4] = 100 * p[4] + 1000  # emission amplitude

    for i in (1, 5):  # abs. and emi. peak positions
        p[i] = 0.05 * p[i] - 0.025 + 8542

    for i in (2, 3, 6, 7):  # Voigt sigmas and gammas
        p[i] = 0.1 * p[i] + 0.1

    p[8] = 300 * p[8] + 2000  # intensity background constant

    # Define each spectrum's classification
    c = [0, 2, 0, 2, 0, 2]
    # Choose single or double component spectrum
    # based on this inside the function `v()`.

    # Generate the spectra
    specs = [v(c[i], w, *p[:, i]) for i in range(6)]

    # Reshape to 2x3 grid
    return np.asarray(specs).reshape((2, 3, len(w)))

raw_data = s()

print('shape of spectral grid:', raw_data.shape)

#%%
# These spectra look as follows,

fig, axes = plt.subplots(2, 3, constrained_layout=True)
for ax, spec in zip(axes.flat, raw_data.reshape((6, raw_data.shape[-1]))):

    plot_spectrum(w, spec, normalised=False, ax=ax)

plt.show()

#%%
# Compute the background intensity
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# `MCALF` does not model a constant background value,
# i.e., the fitting process assumes the intensity
# values far out in the spectrum's wings are zero.
#
# You can however tell `MCALF` what the constant
# background value is, and it will subtract it
# from the spectrum before fitting.
#
# This process was not made automatic as we wanted
# to give the user full control on setting the
# background value.
#
# For these demo data, we shall simply set the background
# to the first intensity value of every spectrum.
# For a real dataset, you may wish to average the
# background value throughout a range of spectral points,
# or even do a moving average throughout time.
# Functions are provided in ``mcalf.utils.smooth`` to
# assist with this process.

backgrounds = raw_data[:, :, 0]

print('shape of background intensity grid:', backgrounds.shape)

#%%
# Define a demo classifier
# ~~~~~~~~~~~~~~~~~~~~~~~~
#
# In this example we will not demonstrate how to create
# a neural network classifier. `MCALF` offers a lot of
# flexibility when it comes to the classifier. This
# demo classifier outlines the basic API that is
# required. By default, the model has a
# :class:`sklearn.neural_network.MLPClassifier`
# object preloaded for use as the classifier.
# The methods :meth:`mcalf.models.ModelBase.train`
# and :meth:`mcalf.models.ModelBase.test` are
# provided by `MCALF` to assist with training the
# neural network. There is also a useful script
# in the `Getting Started` section under
# `User Documentation` in the sidebar for
# semi-automating the process of creating the
# ground truth dataset.
#
# Please see tutorials for packages such
# as scikit-learn for more in-depth advice on
# creating classifiers.
#
# As we only have six spectra with two distinct
# shapes, we can create a very simple classifier
# that classifies spectra based on whether their
# central intensity is greater or smaller than
# the left most intensity.

class DemoClassifier:
    @staticmethod
    def predict(X):
        y = np.zeros(len(X), dtype=int)
        y[X[:, len(X[0]) // 2] > X[:, 0]] = 2
        return y

#%%
# Using ``IBIS8542Model`` with the generated data
# -----------------------------------------------
#
# Initialise the model
# ~~~~~~~~~~~~~~~~~~~~
#
# Everything we have been doing up to this point has been
# creating the demo data and classifier. Now we can actually
# create a model and load in the demo data, although the
# following steps would be identical for a real dataset.

from mcalf.models import IBIS8542Model

# Initialise the model with the wavelength grid
model = IBIS8542Model(original_wavelengths=w)

# Load the spectral shape classifier
model.neural_network = DemoClassifier

# Load the array of spectra and background intensities
model.load_array(raw_data, ['row', 'column', 'wavelength'])
model.load_background(backgrounds, ['row', 'column'])

#%%
# Fit the loaded data
# ~~~~~~~~~~~~~~~~~~~
#
# We have now fully initialised the model. We can now
# call methods to fit the model to the loaded spectra.
# In the following step we fit all the loaded spectra
# and a 1D list of :class:`~mcalf.models.FitResult`
# objects is returned.

fits = model.fit(row=[0, 1], column=range(3))

print(fits)

#%%
# Plot the fitted parameters
# ~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# The individual components of each fit can now be
# plotted separately,

fig, axes = plt.subplots(2, 3, constrained_layout=True)
for ax, fit in zip(axes.flat, fits):

    model.plot_separate(fit, show_legend=False, ax=ax)

plt.show()

#%%
# Extract processed spectra
# ~~~~~~~~~~~~~~~~~~~~~~~~~
#
# As well as fitting spectra, we can call other methods.
# In this step we'll extract the array of loaded spectra.
# However, these spectra have been re-interpolated to a
# new finer wavelength grid. (This grid can be customised
# when initialising the model.)

spectra = model.get_spectra(row=[0, 1], column=range(3))

print('new shape of spectral grid:', spectra.shape)

spectra_1d = spectra.reshape((6, spectra.shape[-1]))

fig, axes = plt.subplots(2, 3, constrained_layout=True)
for ax, spec in zip(axes.flat, spectra_1d):

    plot_spectrum(model.constant_wavelengths, spec,
                  normalised=False, ax=ax)

plt.show()

#%%
# Classify spectra
# ~~~~~~~~~~~~~~~~
#
# We can also classify the loaded spectra and create plots,

classifications = model.classify_spectra(row=[0, 1], column=range(3))

print('classifications are:', classifications)

#%%
# This function plots the spectra grouped
# by classification,

from mcalf.visualisation import plot_classifications

plot_classifications(spectra_1d, classifications.flatten())

#%%
# This function plots a spatial map of
# the classifications on the 2x3 grid,

from mcalf.visualisation import plot_class_map

plot_class_map(classifications)

#%%
# Merge output data
# ~~~~~~~~~~~~~~~~~
#
# The :class:`~mcalf.models.FitResult` objects in the
# ``fits`` 1D list can be merged into a grid.
# Each of these objects can be appended to a
# single :class:`~mcalf.models.FitResults` object.
# This allows for increased data portability,

from mcalf.models import FitResults

# Initialise with the grid shape and num. params.
results = FitResults((2, 3), 8)

# Append each fit to the object
for fit in fits:
    results.append(fit)

#%%
# Now that we have appended all 6 fits,
# we can access the merged data,

print(results.classifications)
print(results.profile)
print(results.success)
print(results.parameters)

#%%
# Calculate velocities
# ~~~~~~~~~~~~~~~~~~~~
#
# And finally, we can calculate Doppler
# velocities for both the quiescent (absorption)
# and active (emission) regimes.
# (The model needs to be given so the stationary
# line core wavelength is available.)

quiescent = results.velocities(model, vtype='quiescent')
active = results.velocities(model, vtype='active')

from mcalf.visualisation import plot_map

fig, axes = plt.subplots(2, constrained_layout=True)
plot_map(quiescent, ax=axes[0])
plot_map(active, ax=axes[1])
plt.show()

#%%
# Export to FITS file
# ~~~~~~~~~~~~~~~~~~~
#
# The :class:`~mcalf.models.FitResults` object can
# also be exported to a FITS file,

results.save('ibis8542model_demo_fit.fits', model)

#%%
# The file has the following structure,

from astropy.io import fits

with fits.open('ibis8542model_demo_fit.fits') as hdul:
    hdul.info()

#%%
# This example outlines the basics of using the
# :class:`mcalf.models.IBIS8542Model` class.
# Typically millions of spectra would need to be
# fitted and processed at the same time.
# For more details on running big jobs and how
# spectra can be fitted in parallel, see
# the `Working with IBIS data` example in the
# gallery.
#
# .. note::
#    Due to limitations with the documentation
#    hosting provider, the examples in the MCALF
#    documentation are computed using a Python
#    based implementation of the Voigt profile
#    instead of the more efficient C version.
#    When this code is run on your own machine
#    it should take much less time than the
#    value reported below.

ext: py

max_stars_repo_name: ConorMacBride/mcalf

max_issues_repo_name: ConorMacBride/mcalf

max_forks_repo_name: ConorMacBride/mcalf

reason: import numpy,from astropy


text:
import numpy as np
import torch
import scipy as sp
from torch import nn
from torch.nn import functional as F
from torch import optim
from tqdm import trange
import math
import matplotlib
from matplotlib import pyplot as plt
import pyaudio
import struct
import moviepy.editor
import os
import wave
import sys

CHUNK = 4096

def chunk2array(raw):
  data = struct.unpack(str(2*CHUNK)+'h', raw)
  return np.array(data).reshape(-1,2)

# CHECKING FOR AUDIOFILE
if not os.path.isfile("parting_ways.wav"):
  moviepy.editor.VideoFileClip("parting_ways.avi").audio.write_audiofile("parting_ways.wav")

# INIT
wf = wave.open("parting_ways.wav", 'rb')
p = pyaudio.PyAudio()
stream = p.open(wf.getframerate(), wf.getnchannels(),p.get_format_from_width(2), output=True)
plt.ion()
fig, ax = plt.subplots()

# START PLOTTING
raw = wf.readframes(CHUNK)
data = np.average(chunk2array(raw), axis=1)
x = np.arange(0, data.size)
line, = ax.plot(x, data)
#ax.set_xlim([xmin,xmax])
ax.set_ylim([-2**15,(2**15)-1])
fig.canvas.draw()
fig.canvas.flush_events()

while True:
  raw = wf.readframes(CHUNK)
  data = np.average(chunk2array(raw), axis=1)
  line.set_ydata(data)
  ax.draw_artist(ax.patch)
  ax.draw_artist(line)
  #fig.canvas.update()
  fig.canvas.flush_events()
  stream.write(raw)

class Net(nn.Module):

  def __init__(self):
    super(Net, self).__init__()
    self.fc1 = nn.Linear(400, 120)  # 6*6 from image dimension
    self.fc2 = nn.Linear(120, 84)
    self.fc3 = nn.Linear(84, 10)

  def forward(self, x):
    # Max pooling over a (2, 2) window
    x = F.max_pool2d(F.relu(self.conv1(x)), (2, 2))
    # If the size is a square you can only specify a single number
    x = F.max_pool2d(F.relu(self.conv2(x)), 2)
    x = x.view(-1, self.num_flat_features(x))
    x = F.relu(self.fc1(x))
    x = F.relu(self.fc2(x))
    x = F.softmax(self.fc3(x), dim=0)
    #x = self.fc3(x)
    return x

base = np.arange(0,100,0.1)[:100]
signal = np.sin(base*6.282)
fourier = np.fft.fft(signal).real[:int(signal.size/2)]
freq = np.fft.fftfreq(base.size,0.05)[:int(signal.size/2)]
print(freq)
plt.plot(freq, fourier, base, signal)
plt.show()

max_stars_repo_name: augustinharter/deeptrip

max_issues_repo_name: augustinharter/deeptrip

max_forks_repo_name: augustinharter/deeptrip

reason: import numpy,import scipy


text:
# Copyright 2021 The Brax Authors.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""Wraps the core environment with some extra statistics for training."""

from typing import Callable, Dict, Tuple

from brax import envs
import flax
import jax
import jax.numpy as jnp

@flax.struct.dataclass
class EnvState:
  """Contains training state for the learner."""
  core: envs.State
  current_episode_metrics: Dict[str, jnp.ndarray]
  completed_episodes_metrics: Dict[str, jnp.ndarray]
  completed_episodes: jnp.ndarray
  completed_episodes_steps: jnp.ndarray

PRNGKey = jnp.ndarray
Action = jnp.ndarray
StepFn = Callable[[EnvState, Action], EnvState]

def wrap(core_env: envs.Env, rng: jnp.ndarray) -> Tuple[EnvState, StepFn]:
  """Returns a wrapped state and step function for training."""
  rng = jax.random.split(rng, core_env.batch_size)

  first_core = core_env.reset(rng)
  first_core.metrics['reward'] = first_core.reward

  def step(state: EnvState, action: Action) -> EnvState:
    core = core_env.step(state.core, action)
    core.metrics['reward'] = core.reward
    # This must be run before test_done in order not to override .steps.
    completed_episodes_steps = state.completed_episodes_steps + jnp.sum(
        core.steps * core.done)
    def test_done(a, b):
      if a is first_core.done or a is first_core.metrics or a is first_core.reward:
        return b
      test_shape = [a.shape[0],] + [1 for _ in range(len(a.shape) - 1)]
      return jnp.where(jnp.reshape(core.done, test_shape), a, b)
    core = jax.tree_multimap(test_done, first_core, core)

    current_episode_metrics = jax.tree_multimap(lambda a, b: a + b,
                                                state.current_episode_metrics,
                                                core.metrics)
    completed_episodes = state.completed_episodes + jnp.sum(core.done)
    completed_episodes_metrics = jax.tree_multimap(
        lambda a, b: a + jnp.sum(b * core.done),
        state.completed_episodes_metrics, current_episode_metrics)
    current_episode_metrics = jax.tree_multimap(
        lambda a, b: a * (1 - core.done) + b * core.done,
        current_episode_metrics, core.metrics)

    return EnvState(
        core=core,
        current_episode_metrics=current_episode_metrics,
        completed_episodes_metrics=completed_episodes_metrics,
        completed_episodes=completed_episodes,
        completed_episodes_steps=completed_episodes_steps)

  first_state = EnvState(
      core=first_core,
      current_episode_metrics=jax.tree_map(jnp.zeros_like, first_core.metrics),
      completed_episodes_metrics=jax.tree_map(
          lambda x: jnp.zeros_like(jnp.sum(x)), first_core.metrics),
      completed_episodes=jnp.zeros(()),
      completed_episodes_steps=jnp.zeros(()))
  return first_state, jax.jit(step)

max_stars_repo_name: benelot/brax

max_issues_repo_name: benelot/brax

max_forks_repo_name: benelot/brax

reason: import jax


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

import gzip
import os
from Bio import SeqIO, pairwise2, motifs
from Bio.Seq import Seq
from datetime import datetime
from scipy import stats
import pandas as pd
import matplotlib.pyplot as plt
import json
from random import randint
import math
import seaborn as sns
import parameters as param

#Reading json from a File
with open('rep_dictionary.json') as json_file:
	rep_dictionary = json.load(json_file)
with open('cont_dictionary.json') as json_file:
	cont_dictionary = json.load(json_file)
with open('total_rep_dictionary.json') as json_file:
	total_rep_dictionary = json.load(json_file)
with open('gold_dictionary.json') as json_file:
	gold_dictionary = json.load(json_file)

# Function to calculate correlation coefficients between two arrays
def corr(x, y, **kwargs):
	# Calculate the values
	coef, p_val = stats.pearsonr(x, y)
	sp, sp_p_val = stats.spearmanr(x, y)
	# Make the label
	label = r'$\rho$ = %s, pearson p-val = %s, sp = %s, spearman p-val = %s.' %(str(round(coef, 2)), str(round(p_val, 2)), str(round(sp, 2)), str(round(sp_p_val, 2)))

	# Add the label to the plot
	ax = plt.gca()
	ax.annotate(label, xy = (0.05, 0.95), size = 20, xycoords = ax.transAxes)

def make_tables(
		uniq_gold_dictionary,
		rep_dictionary,
		cont_dictionary,
		total_rep_dictionary
		):
	mut_dict = {}
	mut_log2_dict = {}
	wt_mut_log2 = {}
	log_2_wt_gaus = {}
	lines_of_text = []
	with open('uniq_gold_dictionary.json') as json_file:
		uniq_gold_dictionary = json.load(json_file)
	for lib in rep_dictionary:
		df_cont_dict = {
				'Barcode': ['control1', 'control2', 'experiment1', 'experiment2'],
				'ROI': ['control1', 'control2', 'experiment1', 'experiment2']
				}

		if uniq_gold_dictionary[lib]['barcodes']:
			if uniq_gold_dictionary[lib]['barcodes-ROIs']:
				df_lib_dict = {'Barcode': [], 'ROI': []}
				for barcode in uniq_gold_dictionary[lib]['barcodes-ROIs']:
					df_lib_dict['Barcode'].append(barcode)
					for ROI in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode]['main ROI']:
						df_lib_dict['ROI'].append(ROI)
			else:
				df_lib_dict = {'Barcode': []}
				for barcode in uniq_gold_dictionary[lib]['barcodes']:
					df_lib_dict['Barcode'].append(barcode)
			unique_barcode_count = len(uniq_gold_dictionary[lib]['barcodes'])

			for rep_type in rep_dictionary[lib]:
				file_count = 0
				if rep_type != 'mapping' and rep_dictionary[lib][rep_type]:
					df_cont_dict[rep_type + ' mean'] = [0, 0, 0, 0]

				if rep_type != 'mapping' and rep_dictionary[lib][rep_type]:
					df_lib_dict[rep_type + ' mean'] = [
						0 for x in range(unique_barcode_count)
						]
					replicate_count = 0
					replicate_list = []

					for file in rep_dictionary[lib][rep_type]:
						file_count += 1
						for replicate in rep_dictionary[lib][rep_type][file]:
							if rep_type != 'mapping':
								replicate_list.append(
									rep_type + ' ' + replicate
									)
							df_dict = {
									'Genuine Barcode': [],
									'Variants of Barcode': [],
									'Count': [],
									'Total Count': []
									}
							df_lib_dict[rep_type + ' ' + replicate] = [
									'' for a in range(unique_barcode_count)
									]

							if rep_type != 'mapping':
								df_cont_dict[rep_type + ' ' + replicate] = [
										'', '', '', ''
										]
							barcode_number = 0
							for cont_barcode in df_cont_dict['Barcode']:
								if rep_type != 'mapping' and cont_barcode in cont_dictionary[lib][rep_type][file][replicate]:
									df_cont_dict[rep_type + ' ' + replicate][barcode_number] = cont_dictionary[lib][rep_type][file][replicate][cont_barcode]
									df_cont_dict[rep_type + ' mean'][barcode_number] += cont_dictionary[lib][rep_type][file][replicate][cont_barcode]
								barcode_number += 1

							for barcode in uniq_gold_dictionary[lib]['barcodes']:
								barcode_index = df_lib_dict['Barcode'].index(
										barcode
										)
								variants_of_barcode = barcode + '\n'
								if barcode in rep_dictionary[lib][rep_type][file][replicate]:
									barcode_count = str(
											rep_dictionary[lib][rep_type][file][replicate][barcode]
											) + '\n'
									total_count = rep_dictionary[lib][rep_type][file][replicate][barcode]
								else:
									barcode_count = '0\n'
									total_count = 0
								for bar_variant in uniq_gold_dictionary[lib]['barcodes'][barcode]:
									variants_of_barcode += bar_variant + '\n'
									if bar_variant in rep_dictionary[lib][rep_type][file][replicate]:
										barcode_count += str(rep_dictionary[lib][rep_type][file][replicate][bar_variant]) + '\n'
										total_count += rep_dictionary[lib][rep_type][file][replicate][bar_variant]
									else:
										barcode_count += '0\n'

								df_dict['Genuine Barcode'].append(barcode)
								df_dict['Variants of Barcode'].append(
										variants_of_barcode
										)
								df_dict['Count'].append(barcode_count)
								df_dict['Total Count'].append(total_count)
								df_lib_dict[rep_type + ' mean'][barcode_index] += total_count
								df_lib_dict[rep_type + ' ' + replicate][barcode_index] = total_count

							df = pd.DataFrame(df_dict, columns=[
									'Genuine Barcode', 'Variants of Barcode',
									'Count', 'Total Count'
									])
							df = df.sort_values(
									by=['Total Count'], ascending=False
									)
							table_path = lib + '_' + rep_type + '_' + replicate + '.xlsx'
							table = df.to_excel(table_path, index = None)
					if replicate_count == 0:
						continue
					df_lib_dict[rep_type + ' mean'] = [
							y / replicate_count for y in df_lib_dict[rep_type + ' mean']
							]
					if rep_type != 'mapping':
						df_cont_dict[rep_type + ' mean'] = [
								d / replicate_count for d in df_cont_dict[rep_type + ' mean']
								]

					scatter_df = pd.DataFrame(data=df_lib_dict)
					#Make joinplot.
					if replicate_count >= 2:
						plt.figure()
						sns.set_style('ticks')
						sns.set_context('poster')
						# Create a pair grid instance
						grid = sns.PairGrid(
								data=scatter_df,
								vars =replicate_list,
								height=10
								)

						# Map the plots to the locations
						grid = grid.map_upper(
								sns.regplot, line_kws={"color": "g"},
								scatter_kws={'s': 0.1}
								)
						grid = grid.map_upper(corr)
						grid = grid.map_lower(
								sns.kdeplot, shade=True
								)
						grid = grid.map_diag(
								sns.kdeplot, linewidth=3,
								shade=True
								)
						plt.savefig(
								'pair_plot_'+ rep_type + lib + '.pdf',
								format='pdf', dpi=1000
								)
						plt.close()

					if total_rep_dictionary[lib][rep_type]:
						df_total_rep_dict = {
									'Genuine Barcode': [],
									'Variants of Barcode': [],
									'Count': [],
									'Total Count': []
									}
						for barcode in uniq_gold_dictionary[lib]['barcodes']:
							var_of_barcode = barcode + '\n'
							barc_count = str(
									total_rep_dictionary[lib][rep_type][barcode] / replicate_count
									) + '\n'
							total_barc_count = total_rep_dictionary[lib][rep_type][barcode] / replicate_count
							for bar_variant in uniq_gold_dictionary[lib]['barcodes'][barcode]:
								var_of_barcode += bar_variant + '\n'
								barc_count += str(
										total_rep_dictionary[lib][rep_type][bar_variant] / replicate_count
										) + '\n'
								total_barc_count += total_rep_dictionary[lib][rep_type][bar_variant] / replicate_count

							df_total_rep_dict['Genuine Barcode'].append(barcode)
							df_total_rep_dict['Variants of Barcode'].append(
									var_of_barcode
									)
							df_total_rep_dict['Count'].append(barc_count)
							df_total_rep_dict['Total Count'].append(
									total_barc_count
									)

						df_total_rep = pd.DataFrame(df_total_rep_dict, columns=[
									'Genuine Barcode', 'Variants of Barcode',
									'Count', 'Total Count'
									])
						df_total_rep = df_total_rep.sort_values(
								by=['Total Count'], ascending=False
								)
						total_rep_path = lib + '_' + rep_type + '.xlsx'
						df_total_rep.to_excel(total_rep_path, index = None)

				else:
					replicate_count = 0
					for file in rep_dictionary[lib][rep_type]:
						file_count += 1
						for replicate in rep_dictionary[lib][rep_type][file]:
							replicate_count += 1
							df_dict = {
									'Genuine Barcode': [],
									'Main ROI': [],
									'Variants of Barcode': [],
									'ROIs' : [],
									'Count': [],
									'Total Count': [],
									'Chimeric molecules, read count': []
									}
							if rep_type != 'mapping':
								df_cont_dict[rep_type + ' ' + replicate] = [
										'', '', '', ''
										]

							barcode_number = 0
							for cont_barcode in df_cont_dict['Barcode']:
								if rep_type != 'mapping' and cont_barcode in cont_dictionary[lib][rep_type][file][replicate]:
									df_cont_dict[rep_type + ' ' + replicate][barcode_number] = cont_dictionary[lib][rep_type][file][replicate][cont_barcode]
									df_cont_dict[rep_type + ' mean'][barcode_number] += cont_dictionary[lib][rep_type][file][replicate][cont_barcode]
								barcode_number += 1

							for barcode in uniq_gold_dictionary[lib]['barcodes-ROIs']:
								total_count = 0
								hybrid_count = 0
								for ROI in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode]['major ROI']:
									main_ROI = ROI
								variants_of_barcode = barcode + '\n'
								ROIs = main_ROI + '\n'
								if barcode in rep_dictionary[lib][rep_type][file][replicate]:
									if main_ROI in rep_dictionary[lib][rep_type][file][replicate][barcode]:
										count = str(rep_dictionary[lib][rep_type][file][replicate][barcode][main_ROI]) + '\n'
										total_count += rep_dictionary[lib][rep_type][file][replicate][barcode][main_ROI]
									else:
										count = '0\n'
								else:
									count = '0\n'
								for ROI_variant in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode][barcode]:
									if ROI_variant != main_ROI:
										variants_of_barcode += '\n'
										ROIs += ROI_variant + '\n'
										if barcode in rep_dictionary[lib][rep_type][file][replicate]:
											if ROI_variant in rep_dictionary[lib][rep_type][file][replicate][barcode]:
												count += str(rep_dictionary[lib][rep_type][file][replicate][barcode][ROI_variant]) + '\n'
												total_count += rep_dictionary[lib][rep_type][file][replicate][barcode][ROI_variant]
												if ROI_variant not in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode]['major ROI'][main_ROI]:
													hybrid_count += rep_dictionary[lib][rep_type][file][replicate][barcode][ROI_variant]
											else:
												count += '0\n'
										else:
											count += '0\n'
								for barc_variant in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode]:
									if len(barc_variant) > 15 and barc_variant != barcode:
										variants_of_barcode += barc_variant
										for mut_variant in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode][barc_variant]:
											variants_of_barcode += '\n'
											ROIs += mut_variant + '\n'
											if barc_variant in rep_dictionary[lib][rep_type][file][replicate]:
												if mut_variant in rep_dictionary[lib][rep_type][file][replicate][barc_variant]:
													count += str(rep_dictionary[lib][rep_type][file][replicate][barc_variant][mut_variant])
													total_count += rep_dictionary[lib][rep_type][file][replicate][barc_variant][mut_variant]
													if mut_variant != main_ROI and mut_variant not in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode]['major ROI'][main_ROI]:
														hybrid_count += rep_dictionary[lib][rep_type][file][replicate][barc_variant][mut_variant]
												else:
													count += '0\n'
											else:
												count += '0\n'
								df_dict['Genuine Barcode'].append(barcode)
								df_dict['Main ROI'].append(
										main_ROI
										)
								df_dict['Variants of Barcode'].append(
										variants_of_barcode
										)
								df_dict['ROIs'].append(ROIs)
								df_dict['Count'].append(count)
								df_dict['Total Count'].append(total_count)
								df_dict['Chimeric molecules, read count'].append(hybrid_count)

							df = pd.DataFrame(df_dict, columns=[
									'Genuine Barcode', 'Main ROI',
									'Variants of Barcode', 'ROIs',
									'Count', 'Total Count', 'Chimeric molecules, read count'
									])
							df['Chimeric molecules, %'] = df['Chimeric molecules, read count'] / df['Total Count'] * 100
							df = df.sort_values(
									by=['Total Count'], ascending=False
									)
							table_path = lib + '_mapping_' + replicate + '.xlsx'
							table = df.to_excel(table_path, index = None)
					if replicate_count == 0:
						continue
					if rep_type != 'mapping':
						df_cont_dict[rep_type + ' mean'] = [
								d / replicate_count for d in df_cont_dict[rep_type + ' mean']
								]

					df_total_map_dict = {
									'Genuine Barcode': [],
									'Main ROI': [],
									'Variants of Barcode': [],
									'ROIs' : [],
									'Count': [],
									'Total Count': [],
									'Chimeric molecules, read count': [],
									'Chimeric molecules, %': []
									}
					for barcode in uniq_gold_dictionary[lib]['barcodes-ROIs']:
						for ROI in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode]['major ROI']:
							main_mut = ROI
						total_map_count = uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode]['total reads'] / replicate_count
						hybr_count = uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode]['hybrid reads'] / replicate_count
						hybrid_percent = hybr_count / total_map_count * 100
						var_of_barc = barcode + '\n'
						mut = main_mut + '\n'
						map_count = str(total_rep_dictionary[lib]['mapping'][barcode][main_mut] / replicate_count) + '\n'
						for mut_var in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode][barcode]:
							if mut_var != main_mut:
								var_of_barc += '\n'
								mut += mut_var + '\n'
								map_count += str(total_rep_dictionary[lib]['mapping'][barcode][mut_var] / replicate_count) + '\n'
						for barc_var in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode]:
							if len(barc_var) > 15 and barc_var != barcode:
								var_of_barcode += barc_var
								for mut_variant in uniq_gold_dictionary[lib]['barcodes-ROIs'][barcode][barc_var]:
									var_of_barc += '\n'
									mut += mut_variant + '\n'
									map_count += str(total_rep_dictionary[lib]['mapping'][barc_var][mut_variant] / replicate_count) + '\n'
						df_total_map_dict['Genuine Barcode'].append(barcode)
						df_total_map_dict['Main ROI'].append(main_mut)
						df_total_map_dict['Variants of Barcode'].append(
								var_of_barc
								)
						df_total_map_dict['ROIs'].append(mut)
						df_total_map_dict['Count'].append(map_count)
						df_total_map_dict['Total Count'].append(
								total_map_count
								)
						df_total_map_dict['Chimeric molecules, read count'].append(hybr_count)
						df_total_map_dict['Chimeric molecules, %'].append(
								hybrid_percent
								)

					df_total_map = pd.DataFrame(df_total_map_dict, columns=[
									'Genuine Barcode', 'Main ROI',
									'Variants of Barcode', 'ROIs',
									'Count', 'Total Count', 'Chimeric molecules, read count',
									'Chimeric molecules, %'
									])
					df_total_map = df_total_map.sort_values(
							by=['Total Count'], ascending=False
							)
					total_map_table_path = lib + '_mapping ' + '.xlsx'
					total_map_table = df_total_map.to_excel(
							total_map_table_path, index = None
							)

			columns_list = [key for key in df_lib_dict]
			df_lib = pd.DataFrame(data=df_lib_dict)
			new_columns_list = ['Barcode']
			pos_start = 1

			if 'ROI' in columns_list:
				new_columns_list.append('ROI')
				pos_start += 1

			new_columns_list += [
					k for k in columns_list if k.startswith('normalization') and not k.endswith('mean')
					]
			new_columns_list += [
					n for n in columns_list if n.startswith('expression') and not n.endswith('mean')
					]
			if 'normalization mean' in columns_list:
				new_columns_list.append('normalization mean')
			if 'expression mean' in columns_list:
				new_columns_list.append('expression mean')
			df_lib = df_lib[new_columns_list]

			cont_new_columns_list = [
					'Barcode', 'ROI'
					] + new_columns_list[pos_start:]

			if 'normalization mean' in new_columns_list and 'expression mean' in new_columns_list:
				df_lib['normalyzed expression'] = df_lib['expression mean'] / df_lib['normalization mean']

				cont_new_columns_list.append('normalyzed expression')
				df_cont_dict['normalyzed expression'] = ['', '', '', '']
				for n in range(4):
					if df_cont_dict['normalization mean'][n] > 0:
						df_cont_dict['normalyzed expression'][n] = df_cont_dict['expression mean'][n] / df_cont_dict['normalization mean'][n]
				if '' not in df_cont_dict['normalyzed expression']:
					cont_new_columns_list.append('normalyzed expression mean')

					df_cont_dict['normalyzed expression mean'] = [
							(df_cont_dict['normalyzed expression'][0] + df_cont_dict['normalyzed expression'][1]) / 2,
							(df_cont_dict['normalyzed expression'][0] + df_cont_dict['normalyzed expression'][1]) / 2,
							(df_cont_dict['normalyzed expression'][2] + df_cont_dict['normalyzed expression'][3]) / 2,
							(df_cont_dict['normalyzed expression'][2] + df_cont_dict['normalyzed expression'][3]) / 2
							]
					if df_cont_dict['normalyzed expression mean'][0] > 0:
						cont_new_columns_list.append('difference')

						if len(param.reference_BCs['forward']) > 0:
							df_lib['normalyzed to reference'] = df_lib['normalyzed expression'] / df_cont_dict['normalyzed expression mean'][0]

						df_cont_dict['difference'] = [
								df_cont_dict['normalyzed expression mean'][2] / df_cont_dict['normalyzed expression mean'][0],
								df_cont_dict['normalyzed expression mean'][2] / df_cont_dict['normalyzed expression mean'][0],
								df_cont_dict['normalyzed expression mean'][2] / df_cont_dict['normalyzed expression mean'][0],
								df_cont_dict['normalyzed expression mean'][2] / df_cont_dict['normalyzed expression mean'][0]
								]

				df_lib = df_lib.sort_values(
						by=['normalyzed expression'], ascending=False
						)

			total_table_path = lib + '.xlsx'
			total_table = df_lib.to_excel(total_table_path, index = None)

			if 'ROI' in new_columns_list:
				#Pairplot for ROIs that have more than two barcodes.
				lib_dict = df_lib.to_dict('index')
				ROI_dictionary = {}
				for index in lib_dict:
					new_ROI = lib_dict[index]['ROI']
					if len(param.reference_BCs['forward']) > 0:
						if new_ROI in ROI_dictionary:
							ROI_dictionary[new_ROI].append(lib_dict[index]['normalyzed to reference'])
						else:
							ROI_dictionary[new_ROI] = [lib_dict[index]['normalyzed to reference']]
					if len(param.reference_BCs['forward']) == 0:
						if new_ROI in ROI_dictionary:
							ROI_dictionary[new_ROI].append(lib_dict[index]['normalyzed expression'])
						else:
							ROI_dictionary[new_ROI] = [lib_dict[index]['normalyzed expression']]

				first_mut_pairplot_list = []
				second_mut_pairplot_list = []

				mut_barcode_count_list = []
				for new_ROI in ROI_dictionary:
					mut_barcode_count_list.append(
							len(ROI_dictionary[new_ROI])
							)
					if len(ROI_dictionary[new_ROI]) >= 2:
						first_random_expr = randint(0, len(ROI_dictionary[new_ROI]) - 1)
						first_mut_pairplot_list.append(ROI_dictionary[new_ROI][first_random_expr])
						second_random_expr = randint(0, len(ROI_dictionary[new_ROI]) - 1)
						while second_random_expr == first_random_expr:
							second_random_expr = randint(0, len(ROI_dictionary[new_ROI]) - 1)
						second_mut_pairplot_list.append(ROI_dictionary[new_ROI][second_random_expr])

				less_11 = 0
				for m in range(1, 11):
					lines_of_text.append(
							'%s ROIs with %s BCs for %s.\n' %(str(mut_barcode_count_list.count(m)), str(m), lib))
					less_11 += mut_barcode_count_list.count(m)
				lines_of_text.append(
						str(len(mut_barcode_count_list) - less_11) + ' ROIs with more than 10 BCs for ' + lib + '.\n')

				plt.figure(figsize=(10,10))
				sns.set_style('ticks')
				sns.set_context('poster')
				sns.displot(
						mut_barcode_count_list, kde=False,
						bins=max(mut_barcode_count_list)
						)
				sns.despine()
				plt.xlabel("Number of BCs for each ROI for " + lib, fontsize=12)
				plt.savefig(
						'BC_per_ROI_' + lib + '.pdf',
						format='pdf', dpi=100
						)
				plt.close()

				plt.figure(figsize=(10,10))
				sns.set_style('ticks')
				sns.set_context('poster')
				sns.displot(
						mut_barcode_count_list, kde=False,
						bins=max(mut_barcode_count_list)
						)
				sns.despine()
				plt.xlim(0, 10)
				plt.xlabel("Number of BCs for each ROI for " + lib, fontsize=12)
				plt.savefig(
						'BC_per_ROI_' + lib + '_limit10' + '.pdf',
						format='pdf', dpi=100
						)
				plt.close()

				pairplot_dict = {
						'first barcode': first_mut_pairplot_list,
						'second barcode': second_mut_pairplot_list
						}
				df_pairplot_dict = pd.DataFrame(
						pairplot_dict,
						columns=['first barcode', 'second barcode']
						)
				plt.figure()
				sns.set_style('ticks')
				sns.set_context('poster')
				# Create a pair grid instance
				grid = sns.PairGrid(
						data=df_pairplot_dict,
						vars = ['first barcode', 'second barcode'],
						height=10
						)

				# Map the plots to the locations
				grid = grid.map_upper(
						sns.regplot, line_kws={"color": "g"},
						scatter_kws={'s': 0.1}
						)
				grid = grid.map_upper(corr)
				grid = grid.map_lower(
						sns.kdeplot, shade=True
						)
				grid = grid.map_diag(
						sns.kdeplot, linewidth=3, shade=True
						)

				plt.savefig(
						'ROI_with_different_BCs_' + lib + '.pdf',
						format='pdf', dpi=1000
						)
				plt.close()

				#Plot mean normalysed expressions for all ROIs.
				if len(param.reference_BCs['forward']) > 0:
					if df_cont_dict['normalyzed expression mean'][0] > 0:
						mut_ser = df_lib.groupby('ROI')['normalyzed to reference'].mean()
					else:
						mut_ser = df_lib.groupby('ROI')['normalyzed expression'].mean()
				if len(param.reference_BCs['forward']) == 0:
					mut_ser = df_lib.groupby('ROI')['normalyzed expression'].mean()
				mut_ser = mut_ser.sort_values(ascending=False)
				mut_table_path = lib + '_ROIs' + '.xlsx'
				mut_table = mut_ser.to_excel(mut_table_path)

				plt.figure(figsize=(10,10))
				sns.set_style('ticks')
				sns.set_context('poster')
				mut_ser_list = mut_ser.tolist()
				sns.displot(
						mut_ser, bins=round(mut_ser_list[0] * 10),
						kde=True
						)
				if len(param.reference_BCs['forward']) > 0:
					if 'normalyzed expression mean' in df_cont_dict:
						plt.axvline(1, 0,1, linewidth=2, color='r')
						plt.axvline(
								df_cont_dict['difference'][0], 0,1, linewidth=2,
								color='g'
								)
				sns.despine()
				plt.savefig(
						'density_plot_' + lib + '.pdf',
						format='pdf', dpi=1000
						)
				if param.wt_ROI[lib] in ROI_dictionary:
					lines_of_text.append(
							'In library %s %s wt ROIs.\n' %(lib, str(len(ROI_dictionary[param.wt_ROI[lib]])))
							)
					sns.displot(
							ROI_dictionary[param.wt_ROI[lib]],
							)
					sns.despine()
					plt.savefig(
						'wt_ROI_normalyzed_expression_' + lib + '.pdf',
						format='pdf', dpi=1000
						)
				plt.close()

				log_2_mut = [math.log2(x) for x in mut_ser_list if x > 0]
				stat, p = stats.shapiro(log_2_mut)
				if p > 0.05:
					lines_of_text.append(
							'For ROIs in lib %s log2 of normalyzed expressions looks Gaussian.\nShapiro p-value=%s.\n' % (lib, str(p))
							)
				else:
					statis, p_value = stats.normaltest(log_2_mut)
					if p_value > 0.05:
						lines_of_text.append(
								'For ROIs in lib %s log2 of normalyzed expressions looks Gaussian.\nShapiro p-value=%s, K^2 p-value=%s.\n' % (lib, str(p), str(p_value))
								)
					else:
						lines_of_text.append(
								'For ROIs in lib %s log2 of normalyzed expressions does not look Gaussian.\nShapiro p-value=%s, K^2 p-value=%s.\n' % (lib, str(p), str(p_value))
								)

				plt.figure(figsize=(10,10))
				sns.set_style('ticks')
				sns.set_context('poster')
				sns.displot(
						log_2_mut, bins=round(log_2_mut[0] * 20),
						kde=True
						)
				if len(param.reference_BCs['forward']) > 0:
					if 'normalyzed expression mean' in df_cont_dict:
						plt.axvline(0, 0,1, linewidth=2, color='r')
						plt.axvline(
								math.log2(df_cont_dict['difference'][0]), 0,1,
								linewidth=2,
								color='g'
								)
				sns.despine()
				plt.savefig(
						'log2_density_plot_' + lib + '.pdf',
						format='pdf', dpi=1000
						)
				if param.wt_ROI[lib] in ROI_dictionary:
					wt_mut_log2[lib] = [math.log2(y) for y in ROI_dictionary[param.wt_ROI[lib]] if y > 0]
					lines_of_text.append('In library %s %s wt ROIs above zero.\n' %(lib, str(len(wt_mut_log2[lib]))))
					statistic, p_val = stats.shapiro(wt_mut_log2[lib])
					if p_val > 0.05:
						lines_of_text.append('For wt in lib %s log2 of normalyzed expressions looks Gaussian.\np-value=%s.\n' % (lib, str(p_val)))
						log_2_wt_gaus[lib] = True
					else:
						lines_of_text.append('For wt in lib %s log2 of normalyzed expressions does not looks Gaussian.\np-value=%s.\n' % (lib, str(p_val)))
						log_2_wt_gaus[lib] = False

					sns.displot(
							wt_mut_log2[lib],
							)
					sns.despine()
					plt.savefig(
						'log2_wt_ROI_normalyzed_expression_' + lib + '.pdf',
						format='pdf', dpi=1000
						)
				plt.close()

				mut_dict[lib] = list(zip(mut_ser.index, mut_ser_list))
				mut_log2_dict[lib] = list(zip(mut_ser.index, log_2_mut))

			df_cont = pd.DataFrame(data=df_cont_dict)
			cont_table_path = lib + '_reference' + '.xlsx'
			if len(param.reference_BCs['forward']) > 0:
				cont_table = df_cont.to_excel(cont_table_path, index = None)

		else:
			for rep_type in rep_dictionary[lib]:
				if len(param.reference_BCs['forward']) > 0:
					df_cont_dict[rep_type + ' mean'] = [0, 0, 0, 0]
					file_count = 0
					for file in rep_dictionary[lib][rep_type]:
						file_count += 1
						for replicate in rep_dictionary[lib][rep_type][file]:
							df_cont_dict[rep_type + ' ' + replicate] = [
										'', '', '', ''
										]
							barcode_number = 0
							for cont_barcode in df_cont_dict['Barcode']:
								if cont_barcode in cont_dictionary[lib][rep_type][file][replicate]:
									df_cont_dict[rep_type + ' ' + replicate][barcode_number] = cont_dictionary[lib][rep_type][file][replicate][cont_barcode]
									df_cont_dict[rep_type + ' mean'][barcode_number] += cont_dictionary[lib][rep_type][file][replicate][cont_barcode]
								barcode_number += 1

						rep_barcodes = {
								lib: {'barcodes': {}, 'barcodes-ROIs': {}}
								}
						rep_barcodes[lib]['barcodes'] = rep_dictionary[lib][rep_type][file][replicate]
						uniq_rep_dictionary = form_uniq_gold_dictionary(
								rep_barcodes
								)
						for barcode in uniq_rep_dictionary[lib]['barcodes']:
							variants_of_barcode = barcode + '\n'
							barcode_count = str(
									rep_dictionary[lib][rep_type][file][replicate][barcode]
									) + '\n'
							total_count = rep_dictionary[lib][rep_type][file][replicate][barcode]
							for bar_variant in uniq_rep_dictionary[lib]['barcodes'][barcode]:
								variants_of_barcode += bar_variant + '\n'
								barcode_count += str(rep_dictionary[lib][rep_type][file][replicate][bar_variant]) + '\n'
								total_count += rep_dictionary[lib][rep_type][file][replicate][bar_variant]
							df_dict['Genuine Barcode'].append(barcode)
							df_dict['Variants of Barcode'].append(
									variants_of_barcode
									)
							df_dict['Count'].append(barcode_count)
							df_dict['Total Count'].append(total_count)
						df = pd.DataFrame(data=df_dict)
						df = df.sort_values(
								by=['Total Count'], ascending=False
								)
						table_path = lib + '_' + rep_type + '_' + replicate + '.xlsx'
						table = df.to_excel(table_path, index = None)

			if len(param.reference_BCs['forward']) > 0:
				df_cont = pd.DataFrame(data=df_cont_dict)
				cont_table_path = lib + '_reference' + '.xlsx'
				cont_table = df_cont.to_excel(cont_table_path, index = None)

	print(lines_of_text)
	with open('run_info.txt', 'a') as info_file:
		info_file.writelines(lines_of_text)

	with open('ROI_dictionary.json', 'w') as f:
		json.dump(mut_dict, f)

	with open('ROI_log2_dictionary.json', 'w') as file:
		json.dump(mut_log2_dict, file)

	with open('wt_log2.json', 'w') as fil:
		json.dump(wt_mut_log2, fil)

	return mut_dict, mut_log2_dict, wt_mut_log2, log_2_wt_gaus

mut_dictionary, mut_log2_dictionary, wt_log2, log_2_wt_gaussian = make_tables(gold_dictionary, rep_dictionary, cont_dictionary, total_rep_dictionary)

def make_graphs(uniq_gold_dictionary):
	for lib in uniq_gold_dictionary:
		if uniq_gold_dictionary[lib]['barcodes']:
			barcode_len_list = [
					len(x) for x in uniq_gold_dictionary[lib]['barcodes']
					]

			plt.figure(figsize=(10,10))
			sns.set_style('ticks')
			sns.set_context('poster')
			sns.displot(
					barcode_len_list, kde=False, bins=4
					)
			sns.despine()
			plt.xlabel("BC length for " + lib, fontsize=14)
			plt.savefig(
					'BC_length_' + lib + '.pdf',
					format='pdf', dpi=1000
					)
			plt.close()

make_graphs(gold_dictionary)

def make_pwm(mut_dict):
	lines_of_text = []
	rounded_mismatch = round(param.ROI_mismatch)
	for lib in mut_dict:
		nucleotide_count = {'A': 0, 'T': 0, 'G': 0, 'C': 0}
		total_count = 0
		for mut_tuple in mut_dict[lib]:
			for nucleotide in mut_tuple[0]:
				nucleotide_count[nucleotide] += 1
				total_count += 1
		lines_of_text.append('ROIs in lib ' + lib + ' contains\n')
		print('ROIs in lib ' + lib + ' contains ')
		for nucl in nucleotide_count:
			nucleotide_count[nucl] = nucleotide_count[nucl] / total_count
			lines_of_text.append(
					str(nucleotide_count[nucl]) + ' of ' + nucl + '.\n'
					)
			print(str(nucleotide_count[nucl]) + ' of ' + nucl)
		for i in range(
				param.ROI_length - rounded_mismatch, param.ROI_length + rounded_mismatch + 1
				):
			instances = [Seq(k[0]) for k in mut_dict[lib] if len(k[0]) == i]
			lines_of_text.append(
					str(len(instances)) + ' ROIs ' + str(i) + ' bp long for ' + lib + '.\n'
					)
			print(str(len(instances)) + ' ROIs ' + str(i) + ' bp long for ' + lib)
			if instances:
				m = motifs.create(instances)
				pwm = m.counts.normalize(pseudocounts=nucleotide_count)
				print('PWM for all ROIs' + str(i) + ' bp long in lib ' + lib)
				print(pwm)
				print(pwm.consensus)
				print(pwm.degenerate_consensus)
				if len(instances) > 20:
					m_high = motifs.create(instances[:round(len(instances) * 0.05)])
					pwm_high = m_high.counts.normalize(pseudocounts=nucleotide_count)
					print('PWM for high expressed ROIs' + str(i) + ' bp long in lib ' + lib)
					print(pwm_high)
					print(pwm_high.consensus)
					print(pwm_high.degenerate_consensus)
					m_low = motifs.create(instances[round(len(instances) * 0.95):])
					pwm_low = m_low.counts.normalize(pseudocounts=nucleotide_count)
					print('PWM for low expressed ROIs' + str(i) + ' bp long in lib ' + lib)
					print(pwm_low)
					print(pwm_low.consensus)
					print(pwm_low.degenerate_consensus)
	with open('run_info.txt', 'a') as info_file:
		info_file.writelines(lines_of_text)

make_pwm(mut_dictionary)

def find_k_mer(mut_dict, len_k_mer):
	k_mer_dict = {}
	k_mer_count = 0
	for mut_tuple in mut_dict:
		ROI = mut_tuple[0]
		for i in range(len(ROI) - len_k_mer + 1):
			k_mer = ROI[i:i + len_k_mer]
			k_mer_count += 1
			if k_mer not in k_mer_dict:
				k_mer_dict[k_mer] = 1
			else:
				k_mer_dict[k_mer] += 1
	for k_mer in k_mer_dict:
		k_mer_dict[k_mer] = k_mer_dict[k_mer] / k_mer_count * 100
	return k_mer_dict

def count_k_mer(ROI_dict):
	lines_of_text = []
	for k_mer_len in range(4, 8):
		for lib in ROI_dict:
			signif_k_mer = {
				'k_mer': [], 'Total count': [], 'High count': [], 'Low count': []
				}
			k_mer_diction = find_k_mer(ROI_dict[lib], k_mer_len)
			lines_of_text.append(
					str(len(k_mer_diction)) + ' k-mers ' + str(k_mer_len) +  ' bp long were found in library ' + lib + '.\n'
					)
			k_mer_high_diction = find_k_mer(
					ROI_dict[lib][:round(len(ROI_dict[lib]) * 0.05)], k_mer_len
					)
			lines_of_text.append(
					str(len(k_mer_high_diction)) + ' high k-mers ' + str(k_mer_len) +  ' bp long were found in library ' + lib + '.\n'
					)
			k_mer_low_diction = find_k_mer(
					ROI_dict[lib][round(len(ROI_dict[lib]) * 0.95):], k_mer_len
					)
			lines_of_text.append(
					str(len(k_mer_low_diction)) + ' low k-mers ' + str(k_mer_len) +  ' bp long were found in library ' + lib + '.\n'
					)
			print(lines_of_text[len(lines_of_text) - 3:])
			for high_k_mer in k_mer_high_diction:
				if k_mer_high_diction[high_k_mer] / (k_mer_len - 1) >= k_mer_diction[high_k_mer]:
					if high_k_mer not in k_mer_low_diction:
						signif_k_mer['k_mer'].append(high_k_mer)
						signif_k_mer['Total count'].append(
								k_mer_diction[high_k_mer]
								)
						signif_k_mer['High count'].append(
								k_mer_high_diction[high_k_mer]
								)
						signif_k_mer['Low count'].append(0)
					elif k_mer_low_diction[high_k_mer] * (k_mer_len - 1) <= k_mer_diction[high_k_mer]:
						signif_k_mer['k_mer'].append(high_k_mer)
						signif_k_mer['Total count'].append(
								k_mer_diction[high_k_mer]
								)
						signif_k_mer['High count'].append(
								k_mer_high_diction[high_k_mer]
								)
						signif_k_mer['Low count'].append(
								k_mer_low_diction[high_k_mer]
								)
			for low_k_mer in k_mer_low_diction:
				if k_mer_low_diction[low_k_mer] / (k_mer_len - 1) >= k_mer_diction[low_k_mer]:
					if low_k_mer not in k_mer_high_diction:
						signif_k_mer['k_mer'].append(low_k_mer)
						signif_k_mer['Total count'].append(
								k_mer_diction[low_k_mer]
								)
						signif_k_mer['High count'].append(0)
						signif_k_mer['Low count'].append(
								k_mer_low_diction[low_k_mer]
								)
					elif k_mer_high_diction[low_k_mer] * (k_mer_len - 1) <= k_mer_diction[low_k_mer]:
						signif_k_mer['k_mer'].append(low_k_mer)
						signif_k_mer['Total count'].append(
								k_mer_diction[low_k_mer]
								)
						signif_k_mer['High count'].append(
								k_mer_high_diction[low_k_mer]
								)
						signif_k_mer['Low count'].append(
								k_mer_low_diction[low_k_mer]
								)
			df = pd.DataFrame(data=signif_k_mer)
			df['High count / Total count'] = df['High count'] / df['Total count']
			df['Low count / Total count'] = df['Low count'] / df['Total count']
			df = df.sort_values(
								by=['High count / Total count'], ascending=False
								)
			table_path = lib + '_k-mer_' + str(k_mer_len) + '.xlsx'
			table = df.to_excel(table_path, index = None)

count_k_mer(mut_dictionary)

def find_log2_k_mer(ROI_lst, len_k_mer):
	k_mer_dict = {}
	for mut_tuple in ROI_lst:
		ROI = mut_tuple[0]
		for i in range(len(ROI) - len_k_mer + 1):
			k_mer = ROI[i:i + len_k_mer]
			if k_mer not in k_mer_dict:
				k_mer_dict[k_mer] = [mut_tuple[1]]
			else:
				k_mer_dict[k_mer].append(mut_tuple[1])
	return k_mer_dict

def log2_k_mer_stat_test(mut_log2_dict, wt_mut_log2, log_2_wt_gaus):
	for k_mer_len in range(4, 7):
		for lib in mut_log2_dict:
			signif_k_mer = {'k_mer': [], 'Mean expression': []}
			k_mer_diction = find_log2_k_mer(mut_log2_dict[lib], k_mer_len)
			for k_mer in k_mer_diction:
				if len(k_mer_diction[k_mer]) >= 20:
					stat, p = stats.shapiro(k_mer_diction[k_mer])
					if p > 0.05 and log_2_wt_gaus[lib]:
						statis, p_value = stats.ttest_ind(
								wt_mut_log2[lib], k_mer_diction[k_mer]
								)
					else:
						statis, p_value = stats.mannwhitneyu(
								wt_mut_log2[lib], k_mer_diction[k_mer]
								)
					if p_value <= 0.05:
						signif_k_mer['k_mer'].append(k_mer)
						k_mer_mean = sum(k_mer_diction[k_mer]) / float(len(k_mer_diction[k_mer]))
						signif_k_mer['Mean expression'].append(k_mer_mean)
			df = pd.DataFrame(data=signif_k_mer)
			df = df.sort_values(by=['Mean expression'], ascending=False)
			table_path = lib + '_stat_k-mer_' + str(k_mer_len) + '.xlsx'
			table = df.to_excel(table_path, index = None)

log2_k_mer_stat_test(mut_log2_dictionary, wt_log2, log_2_wt_gaussian)

reason: from scipy


text:
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras.layers import Flatten, Dense, Conv2D, Dropout, MaxPooling2D
from tensorflow.keras.callbacks import ModelCheckpoint
from tensorflow.keras.models import load_model
from sklearn.model_selection import train_test_split
import cv2

#Initialize Data
data = np.load('/data.npy')
target=np.load('/target.npy')

#Print image with mask on as visual example
plt.imshow(1-data[1100][:, :, 0], cmap='gray');

#Split into train-test sets
X_train, X_test, y_train, y_test=train_test_split(data,target,test_size=0.1,random_state=100)

#Initialize Sequential Model
model = keras.Sequential([
    Conv2D(100, kernel_size=(3,3), input_shape=data.shape[1:]),
    MaxPooling2D(pool_size=(2, 2)),
    Conv2D(100, kernel_size=(3,3)),
    MaxPooling2D(pool_size=(2, 2)),
    Flatten(),
    Dense(75, activation=tf.nn.relu),
    Dense(50, activation=tf.nn.relu),
    Dropout(0.5),
    Dense(2,activation=tf.nn.softmax)
])

#Fit Model
model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])
model.fit(x=X_train,y=y_train, epochs=20)
print(model.evaluate(X_test,y_test))

#Optional: Save model to prevent running
#model.save('/mask_nn.model',overwrite=True,include_optimizer=True,save_format=None,signatures=None,options=None)

#Optional: If model saved, code to load it from Tensorflow rather than running above fitting/compiling section
#model = load_model('/mask_nn.model')

#Test Individual predcitions
image_index = 100
plt.imshow(X_test[image_index][:, :, 0], cmap='gray')
pred = model.predict(X_test[image_index].reshape(1, 100, 100, 1))
print(pred.argmax())

#Test on webcam using OpenCV
clsfr=cv2.CascadeClassifier('/cv2/data/haarcascade_frontalface_default.xml')

source=cv2.VideoCapture(0)
labels_dict={0:'No_Mask',1:'Mask'}
color_dict={0:(0,0,255),1:(0,255,0)}

frame_width = int(source.get(3))
frame_height = int(source.get(4))
out = cv2.VideoWriter('mask_vid.avi',cv2.VideoWriter_fourcc('M','J','P','G'), 10, (frame_width,frame_height))

while(True):
    ret,img=source.read()
    gray=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    face=clsfr.detectMultiScale(gray,1.3,5)

    for x,y,w,h in face:

        face_img=gray[y:y+w,x:x+w]
        normalized=(cv2.resize(face_img,(100,100)))//255.0
        reshaped=np.reshape(normalized,(1,100,100,1))
        result=model.predict(reshaped)

        label=np.argmax(result,axis=1)[0]

        cv2.rectangle(img,(x,y),(x+w,y+h),color_dict[label],2)
        cv2.rectangle(img,(x,y-40),(x+w,y),color_dict[label],-1)
        cv2.putText(
          img, labels_dict[label],
          (x, y-10),
          cv2.FONT_HERSHEY_SIMPLEX,0.8,(255,255,255),2)

    out.write(img)
    cv2.imshow('LIVE',img)

    if cv2.waitKey(1) & 0xFF == ord('q'):
        break

source.release()
out.release()
cv2.destroyAllWindows()

reason: import numpy


text:
import logz
import sys
import matplotlib.pyplot as plt
import cv2
import numpy as np
from kalman import kalman

MOTION_REDUCE = 10500 / 2650
def main(log_filename: str):
    """main."""
    with logz.open(log_filename) as l:
        motion = l["camera_motion"].copy()
        motion_kf = l["camera_motion"].copy()

        motion_kf['x'][:] = kalman(motion['x'], 0.0001, 0.1)[:]

        vrpn = l['vrpn']
        vrpn['x']
        vrpn_dx = logz.dt(vrpn, 'x')
        vrpn_dz = logz.dt(vrpn, 'z')

        # plot for the angles
        plt.subplot(411)
        plt.plot(vrpn['time'], vrpn['x'])

        plt.subplot(412)
        plt.plot(vrpn['time'], vrpn['z'])

        plt.subplot(413)
        plt.plot(motion['time'], motion['x'])

        plt.subplot(414)
        plt.plot(motion['time'], motion['y'])

        plt.show()

if __name__ == '__main__':
    if len(sys.argv) != 2:
        print("usage: %s <run>" % (sys.argv[0],))
        sys.exit(1)
    main(sys.argv[1])

max_stars_repo_name: freddierice/cdrone

max_issues_repo_name: freddierice/cdrone

max_forks_repo_name: freddierice/cdrone


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

import numpy as np
"""
Contains the power function used by MathArray and parsing.py
"""

# This is used in expressions.py's eval_power function also.
def robust_pow(base, exponent):
    """
    Calculates __pow__, and tries other approachs if that doesn't work.

    Usage:
    ======

    >>> robust_pow(5, 2)
    25
    >>> robust_pow(0.5, -1)
    2.0

    If base is negative and power is fractional, complex results are returned:
    >>> almost_j = robust_pow(-1, 0.5)
    >>> np.allclose(almost_j, 1j)
    True
    """
    try:
        return base ** exponent
    except ValueError:
        # In Python 2, (-1)**0.5 raises a ValueError
        # In Python 3, (-1)**0.5 returns the correct value
        return np.lib.scimath.power(base, exponent)


text:
"""
Script that trains TF multitask models on MUV dataset.
"""
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals

import os
import numpy as np
import shutil
import deepchem as dc
from deepchem.molnet import load_muv

np.random.seed(123)

# Load MUV data
muv_tasks, muv_datasets, transformers = load_muv()
train_dataset, valid_dataset, test_dataset = muv_datasets

# Build model
metric = dc.metrics.Metric(
    dc.metrics.roc_auc_score, np.mean, mode="classification")

rate = dc.models.tensorgraph.optimizers.ExponentialDecay(0.001, 0.8, 1000)
model = dc.models.MultiTaskClassifier(
    len(muv_tasks),
    n_features=1024,
    dropouts=[.25],
    learning_rate=rate,
    weight_init_stddevs=[.1],
    batch_size=64,
    verbosity="high")

# Fit trained model
model.fit(train_dataset)

# Evaluate train/test scores
train_scores = model.evaluate(train_dataset, [metric], transformers)
valid_scores = model.evaluate(valid_dataset, [metric], transformers)

print("Train scores")
print(train_scores)

print("Validation scores")
print(valid_scores)


text:
import h5py
import mpi_init
import subfind_data
import cosmology
import tree_build
import numpy as np
from scipy.interpolate import interp1d

def compute_accretion_rates(mpi, path, datasets, data, snap=99, mcut=1.0e12):
    """
    Compute the accretion rate based on different mass estimates
    """

    if not mpi.Rank:
        print(" > Examining: {0}\n > Snap: {1:03d}".format(path, snap), flush=True)

    # First we need to find the halos of interest
    subfind_table = subfind_data.build_table(mpi, sim=path, snap=snap)
    subfind_table.select_halos(mpi, cut=mcut)
    if not mpi.Rank:
        print(" > Found {0:d} halo(s)".format(len(subfind_table.tags)), flush=True)

    # Now rebuild the trees for those halos
    if not mpi.Rank:
        print(" > Building merger tree for halos...", flush=True)
    Mtrees = tree_build.trees(mpi, path, subfind_table, snap)

    Mtrees.build_branches(mpi)

    # Initialize cosmology class instance
    cosmo = cosmology.cosmology(
        subfind_table.hub, subfind_table.omega_m, subfind_table.omega_L
    )

    # Age of Universe at snapshots
    if not mpi.Rank:
        print(" > Computing mass accretion rates", flush=True)
    age_Gyr = cosmo.age(Mtrees.zred)

    # Now compute accretion rates for halos
    tdyn_500c_Gyr = cosmo.t_dynamic_Gyr(Mtrees.zred, delta=500.0, mode="CRIT")
    tdyn_200c_Gyr = cosmo.t_dynamic_Gyr(Mtrees.zred, delta=200.0, mode="CRIT")
    tdyn_200m_Gyr = cosmo.t_dynamic_Gyr(Mtrees.zred, delta=200.0, mode="MEAN")
    tdyn_vir_Gyr = cosmo.t_dynamic_Gyr(Mtrees.zred, mode="VIR")

    # Compute age of Universe one dynamical time ago for each snapshot
    dt_500c_Gyr = age_Gyr - tdyn_500c_Gyr
    dt_200c_Gyr = age_Gyr - tdyn_200c_Gyr
    dt_200m_Gyr = age_Gyr - tdyn_200m_Gyr
    dt_vir_Gyr = age_Gyr - tdyn_vir_Gyr

    # Delta log(a) of a dynamical time for all snapshots
    aexp_int = interp1d(age_Gyr, Mtrees.aexp, fill_value="extrapolate")
    Delta_lgAexp_500c = np.log(Mtrees.aexp) - np.log(aexp_int(dt_500c_Gyr))
    Delta_lgAexp_200c = np.log(Mtrees.aexp) - np.log(aexp_int(dt_200c_Gyr))
    Delta_lgAexp_200m = np.log(Mtrees.aexp) - np.log(aexp_int(dt_200m_Gyr))
    Delta_lgAexp_vir = np.log(Mtrees.aexp) - np.log(aexp_int(dt_vir_Gyr))

    # Now loop over haloes computing Delta log(M) -- with appropriate mass definition
    Delta_lgM500c = np.zeros(Mtrees.M500c.shape, dtype=np.float)
    Delta_lgM200c = np.zeros(Mtrees.M200c.shape, dtype=np.float)
    Delta_lgM200m = np.zeros(Mtrees.M200m.shape, dtype=np.float)
    Delta_lgMvir = np.zeros(Mtrees.Mvir.shape, dtype=np.float)
    for j in range(0, len(Mtrees.index), 1):
        lgM500c_int = interp1d(
            age_Gyr, np.log(Mtrees.M500c[j]), fill_value="extrapolate"
        )
        lgM200c_int = interp1d(
            age_Gyr, np.log(Mtrees.M200c[j]), fill_value="extrapolate"
        )
        lgM200m_int = interp1d(
            age_Gyr, np.log(Mtrees.M200m[j]), fill_value="extrapolate"
        )
        lgMvir_int = interp1d(age_Gyr, np.log(Mtrees.Mvir[j]), fill_value="extrapolate")

        Delta_lgM500c[j] = np.log(Mtrees.M500c[j]) - lgM500c_int(dt_500c_Gyr)
        Delta_lgM200c[j] = np.log(Mtrees.M200c[j]) - lgM200c_int(dt_200c_Gyr)
        Delta_lgM200m[j] = np.log(Mtrees.M200m[j]) - lgM200m_int(dt_200m_Gyr)
        Delta_lgMvir[j] = np.log(Mtrees.Mvir[j]) - lgMvir_int(dt_vir_Gyr)

    # Now compute mass accretion rates
    Macc_500c = Delta_lgM500c / Delta_lgAexp_500c
    Macc_200c = Delta_lgM200c / Delta_lgAexp_200c
    Macc_200m = Delta_lgM200m / Delta_lgAexp_200m
    Macc_vir = Delta_lgMvir / Delta_lgAexp_vir

    # Store and return
    data["Halos"] = Mtrees.halos
    data["Indices"] = Mtrees.index
    data["Redshifts"] = Mtrees.zred
    data["Snapshots"] = Mtrees.snaps
    data["M200mean"] = Mtrees.M200m
    data["M200crit"] = Mtrees.M200c
    data["M500crit"] = Mtrees.M500c
    data["Mvir"] = Mtrees.Mvir
    data["Macc_500c"] = Macc_500c
    data["Macc_200c"] = Macc_200c
    data["Macc_200m"] = Macc_200m
    data["Macc_vir"] = Macc_vir
    del Mtrees
    return

def select_massive_subhalos_with_redshift(mpi, path, data):
    """
    Loop through snapshots reading subfind tables and extracting most massive subhalos
    """

    if not mpi.Rank:
        print(" > Extracting subhalo masses", flush=True)
    masses = data["M200crit"]
    snaps = data["Snapshots"]
    indices = data["Indices"]
    Msubs = np.zeros((len(data["Halos"]), len(snaps), 5), dtype=np.float)
    for j in range(0, len(snaps), 1):
        if not mpi.Rank:
            print("  -{0:2d}".format(snaps[j]), flush=True)
        # First read subfind table to minimum mass
        subfind_table = subfind_data.build_table(mpi, sim=path, snap=snaps[j])
        subfind_table.select_halos(mpi, cut="IDX", max_idx=indices[:, j].max())

        # Now loop over halo indices from the tree trace
        for k in range(0, len(indices[:, j]), 1):
            # Skip halos that are no longer in the tree -- ID == -1
            if indices[k, j] == -1:
                continue
            kdx = np.where(subfind_table.idx == indices[k, j])[0][0]
            mdx = np.argsort(subfind_table.SubMass[kdx])[-5:]
            Msubs[k, j, 0 : len(mdx)] = ((subfind_table.SubMass[kdx])[mdx])[
                ::-1
            ] / 1.989e33
        del subfind_table

    # Store and return
    data["Msubs"] = Msubs
    del masses, snaps, indices, Msubs
    return

def save_data(mpi, path, data):
    """
    Save data to hdf5 file
    """

    tag = {
        "L75n455TNG": "TNG100_L3",
        "L75n910TNG": "TNG100_L2",
        "L75n1820TNG": "TNG100_L1",
        "L205n625TNG": "TNG300_L3",
        "L205n1250TNG": "TNG300_L2",
        "L205n2500TNG": "TNG300_L1",
    }

    # Build output file name
    fname = "output/{0}_tracing.hdf5".format(tag[path.split("/")[-2]])

    # Only Rank zero outputs
    if not mpi.Rank:
        print(" > Saving data to {0}".format(fname), flush=True)
        # Open file
        f = h5py.File(fname, "w")

        # Save some basics
        f.create_dataset("Redshifts", data=data["Redshifts"])
        f.create_dataset("Snapshot_list", data=data["Snapshots"])

        # Loop over halos storing relevant quanities
        for j in range(0, len(data["Halos"]), 1):
            grp = f.create_group(data["Halos"][j])
            f.create_dataset(
                "{0}/FoF_IDs".format(data["Halos"][j]), data=data["Indices"][j]
            )
            f.create_dataset(
                "{0}/M200mean_Msun".format(data["Halos"][j]), data=data["M200mean"][j]
            )
            f.create_dataset(
                "{0}/M200crit_Msun".format(data["Halos"][j]), data=data["M200crit"][j]
            )
            f.create_dataset(
                "{0}/M500crit_Msun".format(data["Halos"][j]), data=data["M500crit"][j]
            )
            f.create_dataset(
                "{0}/Mvir_Msun".format(data["Halos"][j]), data=data["Mvir"][j]
            )
            f.create_dataset(
                "{0}/Macc_500c".format(data["Halos"][j]), data=data["Macc_500c"][j]
            )
            f.create_dataset(
                "{0}/Macc_200c".format(data["Halos"][j]), data=data["Macc_200c"][j]
            )
            f.create_dataset(
                "{0}/Macc_200m".format(data["Halos"][j]), data=data["Macc_200m"][j]
            )
            f.create_dataset(
                "{0}/Macc_vir".format(data["Halos"][j]), data=data["Macc_vir"][j]
            )
            f.create_dataset(
                "{0}/Msubs_Top5_Msun".format(data["Halos"][j]), data=data["Msubs"][j]
            )
        f.close()
    # Wait for save to finish and return
    mpi.comm.Barrier()
    return

if __name__ == "__main__":

    # Simulation of interest
    paths = [
        "/n/hernquistfs3/IllustrisTNG/Runs/L75n455TNG/output",
        "/n/hernquistfs3/IllustrisTNG/Runs/L75n910TNG/output",
        "/n/hernquistfs3/IllustrisTNG/Runs/L75n1820TNG/output",
        "/n/hernquistfs3/IllustrisTNG/Runs/L205n625TNG/output",
        "/n/hernquistfs3/IllustrisTNG/Runs/L205n1250TNG/output",
        "/n/hernquistfs3/IllustrisTNG/Runs/L205n2500TNG/output",
    ]

    # Initialize MPI environment
    mpi = mpi_init.mpi()

    # Datasets required
    datasets = [
        "Group/Group_M_Crit500",
        "Group/Group_M_Crit200",
        "Group/Group_M_Mean200",
        "Group/Group_M_TopHat200",
        "Group/GroupFirstSub",
    ]

    # Loop over simulations
    for x in paths:
        if "L205n" in x:
            cut = 1.0e13
        elif "L75n" in x:
            cut = 1.0e12
        data = {}

        compute_accretion_rates(mpi, x, datasets, data, mcut=cut)

        select_massive_subhalos_with_redshift(mpi, x, data)

        save_data(mpi, x, data)

reason: import numpy,from scipy


text:
import pickle
import numpy as np
import sklearn.datasets
from sklearn.preprocessing import RobustScaler
import pandas as pd
import umap
import matplotlib.pyplot as plt
import sys
import os

def create_umap(X, Y, output_folder=None, layer_str="", save_graph=False, fig=None):

	print(f"Fitting umap")
	reducer = umap.UMAP(random_state=42)

	embedding_train = reducer.fit_transform(X)
	print(f"Done fitting umap")

	if fig is None:
		fig = plt.figure(figsize=(20, 16))

	fig, ax = plt.subplots()
	color_mapping = [int(k) for k in Y]
	scatter = ax.scatter(
		embedding_train[:, 0], embedding_train[:, 1], c=color_mapping, cmap="inferno" , s=10
	)

	legend1 = ax.legend(*scatter.legend_elements(),
			loc="lower left", title="Classes")

	plt.setp(ax, xticks=[], yticks=[])
	plt.suptitle(f"UMAP of RAFT features", fontsize=18)
	ax.set_xlabel(f"UMAP of features")

	plt.show()

	# if save_graph and output_folder is not None:
	# 	save_path = os.path.join(output_folder, f"{layer_str}.png")
	# 	print(f"save_path:{save_path}")
	# 	plt.savefig(save_path, \
	# 		dpi=300, bbox_inches='tight')
	# fig.clf()

path = sys.argv[1]
features = pickle.load(open(path, "rb"))
features = np.array(features)[:100, :]

features = features.reshape(features.shape[0], -1)
create_umap(features, np.arange(len(features)), \
		output_folder='umaps', layer_str='all', save_graph=False)


text:
import glob
import pathlib
import shutil
import subprocess

import femio
import numpy as np

def main():
    original_data_paths = [
        pathlib.Path(cnt_file).parent for cnt_file
        in glob.glob(
            'tests/data/simple/rotated_thermal/raw/**/thermal.cnt',
            recursive=True)]

    for original_data_path in original_data_paths:
        transform_data(original_data_path, additional_trial=2)
    return

def transform_data(data_path, *, additional_trial=2):
    fem_data = femio.FEMData.read_directory('fistr', data_path, read_npy=False)

    orthogonal_matrix = np.array([[1, 0, 0], [0, 0, 1], [0, -1, 0]])
    output_directory = data_path.parent / (data_path.name + 'rotated_x')
    process(fem_data, orthogonal_matrix, output_directory)

    for i in range(additional_trial):
        orthogonal_matrix = generate_rotation_matrix()
        output_directory = data_path.parent / (data_path.name + f"rotated_{i}")
        process(fem_data, orthogonal_matrix, output_directory)
    return

def process(fem_data, rotation_matrix, output_directory):
    _rotate_data(fem_data, rotation_matrix, output_directory)
    sp = subprocess.run(
        f"cd {output_directory} && fistr1", shell=True, check=True)
    print(sp)
    validate_results(fem_data, output_directory, rotation_matrix)
    np.savetxt(output_directory / 'rotation_matrix.txt', rotation_matrix)
    return

def generate_rotation_matrix():
    vec1 = normalize(np.random.rand(3)*2 - 1)
    vec2 = normalize(np.random.rand(3)*2 - 1)
    vec3 = normalize(np.cross(vec1, vec2))
    vec2 = np.cross(vec3, vec1)
    return np.array([vec1, vec2, vec3])

def normalize(x):
    return x / np.linalg.norm(x)

def validate_results(original_fem_data, output_directory, rotation_matrix):
    calculated_fem_data = femio.FEMData.read_directory(
        'fistr', output_directory, read_npy=False)
    rotated_original_strain = rotate_tensor_array(
        original_fem_data,
        original_fem_data.elemental_data.get_attribute_data(
            'ElementalSTRAIN'), rotation_matrix)
    calculated_strain = calculated_fem_data.elemental_data.get_attribute_data(
        'ElementalSTRAIN')
    mean_rmse = np.mean((
        rotated_original_strain - calculated_strain)**2)**.5
    ref = np.mean(rotated_original_strain**2)**.5
    print('========================')
    print(f"mean error: {mean_rmse}")
    print(f"relative mean error: {mean_rmse / ref * 100}")
    print('========================')
    with open(output_directory / 'log.txt', 'w') as f:
        f.write(f"mean error: {mean_rmse}\n")
        f.write(f"relative mean error: {mean_rmse / ref * 100}\n")
    if mean_rmse / ref * 100 > 1e-5:
        raise ValueError('Error too big')
    return

def rotate_tensor_array(fem_data, tensor_array, rotation_matrix):
    symmetric_mat = fem_data.convert_array2symmetric_matrix(
        tensor_array, from_engineering=True)
    rotated_mat = np.array([
        rotation_matrix @ l @ rotation_matrix.T for l in symmetric_mat])
    rotated_array = fem_data.convert_symmetric_matrix2array(
        rotated_mat, to_engineering=True)
    return rotated_array

def _rotate_data(fem_data, rotation_matrix, output_directory):
    shutil.rmtree(output_directory, ignore_errors=True)
    new_fem_data = femio.FEMData(
        fem_data.nodes, fem_data.elements)

    # Nodal data
    original_node = fem_data.nodal_data.get_attribute_data('NODE')
    rotated_nodes = np.array([
        rotation_matrix @ n for n in original_node])
    new_fem_data.nodes.data = rotated_nodes
    original_t_init = fem_data.nodal_data.get_attribute_data(
        'INITIAL_TEMPERATURE')
    original_t_cnt = fem_data.nodal_data.get_attribute_data('CNT_TEMPERATURE')
    nodal_data_dict = {
        'NODE': rotated_nodes,
        'INITIAL_TEMPERATURE': original_t_init,
        'CNT_TEMPERATURE': original_t_cnt}
    new_fem_data.nodal_data.update_data(
        new_fem_data.nodes.ids, nodal_data_dict, allow_overwrite=True)

    # Material data
    original_poisson_ratio = np.mean(
        fem_data.elemental_data.get_attribute_data('Poisson_ratio'),
        axis=0, keepdims=True)
    original_young_modulus = np.mean(
        fem_data.elemental_data.get_attribute_data('Young_modulus'),
        axis=0, keepdims=True)
    original_ltec = np.mean(
        fem_data.elemental_data.get_attribute_data(
            'linear_thermal_expansion_coefficient_full'),
        axis=0, keepdims=True)
    rotated_lte_array = rotate_tensor_array(
        new_fem_data, original_ltec, rotation_matrix)
    material_data_dict = {
        'Poisson_ratio': original_poisson_ratio,
        'Young_modulus': original_young_modulus,
        'linear_thermal_expansion_coefficient_full': rotated_lte_array}
    new_fem_data.materials.update_data('MAT_ALL', material_data_dict)

    # Elemental data
    n_element = len(new_fem_data.elements.ids)
    elemental_data_dict = {
        'Poisson_ratio': original_poisson_ratio * np.ones((n_element, 1)),
        'Young_modulus': original_young_modulus * np.ones((n_element, 1)),
        'linear_thermal_expansion_coefficient_full':
        rotated_lte_array * np.ones((n_element, 6))}
    new_fem_data.elemental_data.update_data(
        new_fem_data.elements.ids, elemental_data_dict)

    # Other info
    new_fem_data.settings = {
        'solution_type': 'STATIC',
        'output_res': 'NSTRAIN,ON\nNSTRESS,ON\n',
        'output_vis': 'NSTRAIN,ON\nNSTRESS,ON\n'}
    new_fem_data.element_groups = {'ALL': fem_data.elements.ids}
    new_fem_data.sections.update_data(
        'MAT_ALL', {'TYPE': 'SOLID', 'EGRP': 'ALL'})
    new_fem_data.constraints['spring'] = fem_data.constraints['spring']
    new_fem_data.material_overwritten = False

    new_fem_data.write('fistr', output_directory / 'mesh')
    return

if __name__ == '__main__':
    main()


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Nov 27 16:58:52 2020

This source code is adapted from SentEval dataset https://github.com/facebookresearch/SentEval.
The code allows you to evaluate Sentence Transformers models in different tasks concatenating
the sentence embeddings from different multilingual models and applying a PCA:

    - distiluse-base-multilingual-cased
    - xlm-r-bert-base-nli-stsb-mean-tokens
    - xlm-r-distilroberta-base-paraphrase-v1
    - LaBSE
    - distilbert-multilingual-nli-stsb-quora-ranking

If you want to use it, you have to clone this SentEval repository and download the
data as it is explained in the examples:
    $  git clone https://github.com/Huertas97/SentEval.git
    $  cd SentEval/data/downstream
    $  bash ./get_transfer_data.bash

For more information about how to use this code read help information:
    $  python ensemble_dim_reduc.py --help

@author: Álvaro Huertas García
"""

from __future__ import absolute_import, division, unicode_literals

import sys
import os
import torch
import logging
import numpy as np
from sentence_transformers import SentenceTransformer
import sys
import pickle
import pandas as pd
from optparse import OptionParser
# Process command-line options
parser = OptionParser(add_help_option=False)

"""
Command line options
"""

# General options
parser.add_option('-d', '--data', type="string", help='Number of sentences to retrieve')
parser.add_option('-o', '--output', type="string", help='Output name')
parser.add_option('-m', '--models', type = "string", help='Languages from which we extract sentences')
parser.add_option('-h', '--help', action='store_true', help='Show this help message and exit.')
(options, args) = parser.parse_args()

def print_usage():
    print("""
Usage:

    python ensembe_dim_reduc.py [options]

Options:
    -d, --data              Path to dataframe with model names and PCAs
    -m, --models            Name models used for SentEval
    -o, --output            Output name
    -h, --help              Help information

Example. Evaluate distiluse-base-multilingual-cased and xlm-r-bert-base-nli-stsb-mean-tokens ensemble on SentEval tasks :
    python ensemble_dim_reduc.py --data df_multi_selected_99.pkl --models distiluse-base-multilingual-cased,xlm-r-bert-base-nli-stsb-mean-tokens""")
    sys.exit()

if not options.data or not options.models:
    print_usage()

if not options.output:
    options.output = options.models

if options.help:
    print_usage()

"""
Class to create a ensemble object with PCA for each model
"""

class ensemble_stransformer:

  def __init__(self, model_names, pca_list):

    # Set models and PCAs as attributes
    for i, (model_arg, pca_arg) in enumerate(zip(model_names, pca_list)):
      sentence_model = SentenceTransformer(model_arg)
      new_model_att = {"model_"+str(i): sentence_model,
                        "model_" + str(i) + "_pca": pca_arg}
      self.__dict__.update(new_model_att)

  def PCA(self, name, embedding):
    for att in dir(self):
      if name + "_pca" == att:
        emb_dim_red = getattr(self, att).transform(embedding)
        return emb_dim_red

  def encode(self, sentences):
    embeddings = []
    n_models = 0
    for att in dir(self):
      if "model_" in att and "pca" not in att:
        emb = getattr(self, att).encode(sentences)
        emb_dim_red = self.PCA(att, emb)
        embeddings.append(emb_dim_red)
        n_models += 1

    else:
      embeddings_concat = np.concatenate( embeddings, axis = 1)
      return embeddings_concat

# Set PATHs
PATH_TO_SENTEVAL = '../'
PATH_TO_DATA = '../data'

# import SentEval
sys.path.insert(0, PATH_TO_SENTEVAL)
import senteval

"""
Prepare and batcher function required by SentEval
"""

# SentEval prepare and batcher
def prepare(params, samples):
    return

def batcher(params, batch):
     sentences = [' '.join(s) for s in batch]
     embeddings = params.ensemble.encode(sentences)
     return embeddings

"""
Evaluation of trained model on Transfer Tasks (SentEval)
"""

# define senteval params
params_senteval = {'task_path': PATH_TO_DATA, 'usepytorch': True, 'kfold': 5}
params_senteval['classifier'] = {'nhid': 0, 'optim': 'rmsprop', 'batch_size': 128,
                                 'tenacity': 3, 'epoch_size': 2}

# Set up logger
logging.basicConfig(format='%(asctime)s : %(message)s', level=logging.DEBUG)

# Extract PCA fitted for the multilingual models
df_pca_data = os.path.join(PATH_TO_DATA, "PCA", options.data)
df_ensemble_pca = pd.read_pickle(df_pca_data)
models_names = df_ensemble_pca.Modelo.to_list()
models_pca = df_ensemble_pca.PCA.to_list()

# Selecting PCA from models introduced by command line
model_names = options.models.split(",")
pca_list = []
for name in model_names:
  pca = df_ensemble_pca[df_ensemble_pca["Modelo"] == name].PCA.values[0]
  pca_list.append(pca)

# Creating the models with PCA that will be used to obtain embeddings by SentEval
params_senteval["ensemble"] = ensemble_stransformer( model_names, pca_list)

# SentEval Engine
if __name__ == "__main__":

    se = senteval.engine.SE(params_senteval, batcher, prepare)
    transfer_tasks = [
                      'CR',
                      'MR', 'MPQA', 'SUBJ', 'SST2', 'SST5', 'TREC', 'MRPC',
                      # 'SNLI',
                      'SICKEntailment', 'SICKRelatedness', 'STSBenchmark',
                      #  'ImageCaptionRetrieval',
                      'STS12', 'STS13', 'STS14', 'STS15', 'STS16',
                      'Length', 'WordContent', 'Depth', 'TopConstituents','BigramShift', 'Tense',
                      'SubjNumber', 'ObjNumber', 'OddManOut', 'CoordinationInversion'
                      ]
    results = se.eval(transfer_tasks)

    # create result folder
    logging.debug("****** Saving results ******")
    result_path = "../result"
    if not os.path.exists(result_path):
      os.mkdir(result_path)

    save_path = os.path.join(result_path, options.output+"_result.pkl")
    pickle.dump( results, open( save_path , "wb" ) )
    logging.debug("Results saved in {}".format(save_path))
    print(results)