<commit_before>#!/usr/bin/env python

from cpyMSpec import centroidize
from pyimzml.ImzMLParser import ImzMLParser
from pyimzml.ImzMLWriter import ImzMLWriter
import numpy as np
import cPickle

import argparse

parser = argparse.ArgumentParser(description="add noise to a clean dataset")
parser.add_argument('real', type=str, help="original centroided imzML")
parser.add_argument('simclean', type=str, help="input file produced by simulateClean.py")
parser.add_argument('nmf', type=str, help="factorization file produced by NNMF.py")
parser.add_argument('layers', type=str, help="layers file produces by assignMolecules.py")
parser.add_argument('output', type=str, help="output filename (centroided .imzML)")

args = parser.parse_args()

class NoiseGenerator(object):
    def __init__(self, nmf_fn, layers_fn, imzml_fn):
        self._imzml = ImzMLParser(imzml_fn)
        with np.load(nmf_fn) as data:
            nx, ny = data['shape']
            self._W = data['W'].reshape((nx, ny, -1))
            self._H = data['H']
            self._mz_axis = data['mz_axis']
        self._coords = {}
        for i, coords in enumerate(self._imzml.coordinates):
            self._coords[(coords[0], coords[1])] = i
        self._mz_bins = []
        for mz, ppm in self._mz_axis:
            self._mz_bins.append(mz * (1.0 + 1e-6 * ppm))

        # self._removeAssignedBins(layers_fn)

    def _removeAssignedBins(self, layers_fn):
        # buggy at the moment
        with open(layers_fn, 'rb') as f:
            layers = cPickle.load(f)
        for i in layers['layers_list']:
            assigned = layers['layers_list'][i]['assigned_mz_bins']
            assigned = assigned[assigned < self._H[i].shape[0]]
            print "#assigned bins in component #{}: {}".format(i + 1, len(assigned))
            h = np.zeros_like(self._H[i])
            h[assigned] = self._H[i][assigned]
            self._H[i] = h

    def _getRealSpectrum(self, x, y):
        return self._imzml.getspectrum(self._coords[(x, y)])

    def generateNoise(self, x, y):
        real_spectrum = self._getRealSpectrum(x, y)
        real_mzs, real_intensities = map(np.array, real_spectrum)

        min_mz, max_mz = self._mz_bins[0], self._mz_bins[-1]
        inside_range = (real_mzs >= min_mz) & (real_mzs <= max_mz)
        real_mzs = real_mzs[inside_range]
        real_intensities = real_intensities[inside_range]

        bins = np.digitize(real_mzs, self._mz_bins)
        n_bins = len(self._mz_bins)
        binned_real_intensities = np.bincount(bins, real_intensities, n_bins)
        binned_approx_intensities = self._W[x, y, :].dot(self._H)
        noise = np.abs(binned_real_intensities - binned_approx_intensities)
        # FIXME: avoid duplicating noise
        noise_intensities = noise[bins]
        noise_mzs = np.array(real_mzs)
        nnz = noise_intensities > min(real_intensities) / 2
        return noise_mzs[nnz], noise_intensities[nnz]

    def addNoise(self, profile_spectrum, coords):
        spec = map(np.array, profile_spectrum)
        p = centroidize(*spec)
        mzs = np.array(p.masses)
        intensities = np.array(p.abundances) * spec[1].max()

        noise_mzs, noise_intensities = self.generateNoise(*coords)
        mzs = np.concatenate([mzs, noise_mzs])
        intensities = np.concatenate([intensities, noise_intensities])<commit_msg>bugfix for empty profile spectra<commit_after>#!/usr/bin/env python

from cpyMSpec import centroidize
from pyimzml.ImzMLParser import ImzMLParser
from pyimzml.ImzMLWriter import ImzMLWriter
import numpy as np
import cPickle

import argparse

parser = argparse.ArgumentParser(description="add noise to a clean dataset")
parser.add_argument('real', type=str, help="original centroided imzML")
parser.add_argument('simclean', type=str, help="input file produced by simulateClean.py")
parser.add_argument('nmf', type=str, help="factorization file produced by NNMF.py")
parser.add_argument('layers', type=str, help="layers file produces by assignMolecules.py")
parser.add_argument('output', type=str, help="output filename (centroided .imzML)")

args = parser.parse_args()

class NoiseGenerator(object):
    def __init__(self, nmf_fn, layers_fn, imzml_fn):
        self._imzml = ImzMLParser(imzml_fn)
        with np.load(nmf_fn) as data:
            nx, ny = data['shape']
            self._W = data['W'].reshape((nx, ny, -1))
            self._H = data['H']
            self._mz_axis = data['mz_axis']
        self._coords = {}
        for i, coords in enumerate(self._imzml.coordinates):
            self._coords[(coords[0], coords[1])] = i
        self._mz_bins = []
        for mz, ppm in self._mz_axis:
            self._mz_bins.append(mz * (1.0 + 1e-6 * ppm))

        # self._removeAssignedBins(layers_fn)

    def _removeAssignedBins(self, layers_fn):
        # buggy at the moment
        with open(layers_fn, 'rb') as f:
            layers = cPickle.load(f)
        for i in layers['layers_list']:
            assigned = layers['layers_list'][i]['assigned_mz_bins']
            assigned = assigned[assigned < self._H[i].shape[0]]
            print "#assigned bins in component #{}: {}".format(i + 1, len(assigned))
            h = np.zeros_like(self._H[i])
            h[assigned] = self._H[i][assigned]
            self._H[i] = h

    def _getRealSpectrum(self, x, y):
        return self._imzml.getspectrum(self._coords[(x, y)])

    def generateNoise(self, x, y):
        real_spectrum = self._getRealSpectrum(x, y)
        real_mzs, real_intensities = map(np.array, real_spectrum)

        min_mz, max_mz = self._mz_bins[0], self._mz_bins[-1]
        inside_range = (real_mzs >= min_mz) & (real_mzs <= max_mz)
        real_mzs = real_mzs[inside_range]
        real_intensities = real_intensities[inside_range]

        bins = np.digitize(real_mzs, self._mz_bins)
        n_bins = len(self._mz_bins)
        binned_real_intensities = np.bincount(bins, real_intensities, n_bins)
        binned_approx_intensities = self._W[x, y, :].dot(self._H)
        noise = np.abs(binned_real_intensities - binned_approx_intensities)
        # FIXME: avoid duplicating noise
        noise_intensities = noise[bins]
        noise_mzs = np.array(real_mzs)
        nnz = noise_intensities > min(real_intensities) / 2
        return noise_mzs[nnz], noise_intensities[nnz]

    def addNoise(self, profile_spectrum, coords):
        spec = map(np.array, profile_spectrum)
        p = centroidize(*spec)
        mzs = np.array(p.masses)
        mult = spec[1].max() if len(spec[1]) > 0 else 1
        intensities = np.array(p.abundances) * mult

        noise_mzs, noise_intensities = self.generateNoise(*coords)
        mzs = np.concatenate([mzs, noise_mzs])
        intensities = np.concatenate([intensities, noise_intensities])
{'fraction_non_alphanumeric': 0.12141802067946825, 'fraction_numerical': 0.004283604135893649, 'mean_word_length': 3.1036363636363635, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2163631', 'n_tokens_mistral': 2414, 'n_tokens_neox': 2369, 'n_words': 468}