Analyze Kmer spectrum, calculations derived from
        allpathslg/src/kmers/KmerSpectra.cc
348987
def analyze(self, ploidy=2, K=23, covmax=1000000):
        
        from math import sqrt

        data = self.data
        kf_ceil = max(K for (K, c) in data)
        if kf_ceil > covmax:
            exceeds = sum(1 for (K, c) in data if K > covmax)
            logging.debug("A total of {0} distinct K-mers appear > "
                          "{1} times. Ignored ...".format(exceeds, covmax))
            kf_ceil = covmax

        nkf = kf_ceil + 1
        a = [0] * nkf
        for kf, c in data:
            if kf > kf_ceil:
                continue
            a[kf] = c

        ndk = a  # number of distinct kmers
        nk = [k * c for k, c in enumerate(a)]  # number of kmers
        cndk = [0] * nkf  # cumulative number of distinct kmers
        cnk = [0] * nkf  # cumulative number of kmers
        for kf in xrange(1, nkf):
            cndk[kf] = cndk[kf - 1] + .5 * (ndk[kf - 1] + ndk[kf])
            cnk [kf] = cnk [kf - 1] + .5 * (nk [kf - 1] + nk [kf])

        # Separate kmer spectrum in 5 regions based on the kf
        # 1        ... kf_min1    : bad kmers with low frequency
        # kf_min1  ... kf_min2    : good kmers CN = 1/2 (SNPs)
        # kf_min2  ... kf_min3    : good kmers CN = 1
        # kf_min3  ... kf_hi      : good kmers CN > 1 (repetitive)
        # kf_hi    ... inf        : bad kmers with high frequency

        # min1: find first minimum
        _kf_min1 = 10
        while (_kf_min1 - 1 >= 2 and nk[_kf_min1 - 1] < nk[_kf_min1]):
            _kf_min1 -= 1
        while (_kf_min1 <= kf_ceil and nk [_kf_min1 + 1] < nk[_kf_min1]):
            _kf_min1 += 1

        # max2: find absolute maximum mx2 above first minimum min1
        _kf_max2 = _kf_min1
        for kf in xrange(_kf_min1 + 1, int(0.8 * kf_ceil)):
            if nk[kf] > nk[_kf_max2]:
                _kf_max2 = kf

        # max2: resetting max2 for cases of very high polymorphism
        if ploidy == 2:
            ndk_half = ndk[_kf_max2 / 2]
            ndk_double = ndk[_kf_max2 * 2]
            if ndk_double > ndk_half:
                _kf_max2 *= 2

        # max1: SNPs local maximum max1 as half global maximum max2
        _kf_max1 = _kf_max2 / 2

        # min2: SNPs local minimum min2 between max1 and max2
        _kf_min2 = _kf_max1 * (2 * ndk[_kf_max1] + ndk[_kf_max2]) / \
                    (ndk[_kf_max1] + ndk[_kf_max2])

        # min1: refine between min1 and max2/2
        for kf in xrange(_kf_min1 + 1, _kf_max1):
            if (nk[kf] < nk[_kf_min1]):
                _kf_min1 = kf

        # min3: not a minimum, really. upper edge of main peak
        _kf_min3 = _kf_max2 * 3 / 2

        print("kfs:", _kf_min1, _kf_max1, \
                        _kf_min2, _kf_max2, _kf_min3, file=sys.stderr)
        self.min1 = _kf_min1
        self.max1 = _kf_max1
        self.min2 = _kf_min2
        self.max2 = _kf_max2
        self.min3 = _kf_min3

        # Define maximum kf above which we neglect data
        _kf_hi = _kf_max2 * sqrt(4 * ndk[2 * _kf_max2] * _kf_max2) \
                    if 2 * _kf_max2 < len(ndk) else \
                 _kf_max2 * sqrt(4 * ndk[len(ndk) - 1] * _kf_max2)
        _kf_hi = int(_kf_hi)

        if _kf_hi > kf_ceil:
            _kf_hi = kf_ceil

        _nk_total = cnk[len(cnk) - 1]
        _nk_bad_low_kf = cnk[_kf_min1]
        _nk_good_uniq = cnk[_kf_min3] - cnk[_kf_min2]
        _nk_bad_high_kf = _nk_total - cnk[_kf_hi]
        _ndk_good_snp = cndk[_kf_min2] - cndk[_kf_min1]
        _ndk_good_uniq = cndk[_kf_min3] - cndk[_kf_min2]

        # kmer coverage C_k
        _kf_ave_uniq = _nk_good_uniq * 1. / _ndk_good_uniq
        _genome_size = (_nk_total - _nk_bad_low_kf - _nk_bad_high_kf) / \
                        _kf_ave_uniq
        _genome_size_unique = _ndk_good_uniq + _ndk_good_snp / 2
        _genome_size_repetitive = _genome_size - _genome_size_unique
        _coverage = _nk_total / _genome_size if _genome_size else 0

        # SNP rate estimation, assumes uniform distribution of SNPs over the
        # genome and accounts for the reduction in SNP kmer counts when
        # polymorphism is very high
        if ploidy == 2:
            _d_SNP = 1. / (1. - (1. - .5 * _ndk_good_snp / _genome_size) ** (1. / K)) \
                     if _ndk_good_snp > 0 else 1000000

        G = int(_genome_size)
        G1 = int(_genome_size_unique)
        GR = int(_genome_size_repetitive)
        coverage = int(_coverage)

        m = "Kmer (K={0}) Spectrum Analysis\n".format(K)
        m += "Genome size estimate = {0}\n".format(thousands(G))
        m += "Genome size estimate CN = 1 = {0} ({1})\n".format(thousands(G1),
                percentage(G1, G))
        m += "Genome size estimate CN > 1 = {0} ({1})\n".format(thousands(GR),
                percentage(GR, G))
        m += "Coverage estimate: {0} x\n".format(coverage)
        self.repetitive = "Repeats: {0} percent".format(GR * 100 / G)

        if ploidy == 2:
            d_SNP = int(_d_SNP)
            self.snprate = "SNP rate ~= 1/{0}".format(d_SNP)
        else:
            self.snprate = "SNP rate not computed (Ploidy = {0})".format(ploidy)
        m += self.snprate + '\n'

        self.genomesize = int(round(self.totalKmers * 1. / self.max2))

        print(m, file=sys.stderr)
1185273
def spectral_centroid_and_spread(data, fs):
    """"""
    data = np.mean(data, axis=1)

    nFFT = len(data) // 2
    X = FFT(data, nFFT)

    ind = (np.arange(1, len(X) + 1)) * (fs / (2.0 * len(X)))

    Xt = X.copy()
    Xt = Xt / Xt.max()
    NUM = np.sum(ind * Xt)
    DEN = np.sum(Xt) + EPSILON

    # Centroid:
    C = (NUM / DEN)

    # Spread:
    S = np.sqrt(np.sum(((ind - C) ** 2) * Xt) / DEN)

    # Normalize:
    C = C / (fs / 2.0)
    S = S / (fs / 2.0)

    return (C, S)
883695
def _extract_peaks(specgram, neighborhood, threshold):
    
    kernel = np.ones(shape=neighborhood)
    local_averages = convolve(specgram, kernel / kernel.sum(), mode="constant", cval=0)

    # suppress all points below the floor value
    floor = (1 + threshold) * local_averages
    candidates = np.where(specgram > floor, specgram, 0)

    # grayscale dilation is equivalent to non-maximal suppression
    local_maximums = grey_dilation(candidates, footprint=kernel)
    peak_coords = np.argwhere(specgram == local_maximums)
    peaks = zip(peak_coords[:, 0], peak_coords[:, 1])

    return peaks
863096
def _reads_in_peaks(bam_file, peaks_file, sample):
    """"""
    if not peaks_file:
        return {}
    rip = number_of_mapped_reads(sample, bam_file, bed_file = peaks_file)
    return {"metrics": {"RiP": rip}}
587371
def _multitaper_cross_spectrum(self, clm, slm, k, convention='power',
                                   unit='per_l', lmax=None, taper_wt=None):
        
        if lmax is None:
            lmax = min(clm.lmax, slm.lmax)

        sh1 = clm.to_array(normalization='4pi', csphase=1, lmax=lmax)
        sh2 = slm.to_array(normalization='4pi', csphase=1, lmax=lmax)

        if taper_wt is None:
            mtse, sd = _shtools.SHMultiTaperMaskCSE(sh1, sh2, self.tapers,
                                                    lmax=lmax, k=k)
        else:
            mtse, sd = _shtools.SHMultiTaperMaskCSE(sh1, sh2, self.tapers,
                                                    lmax=lmax, k=k,
                                                    taper_wt=taper_wt)

        if (unit == 'per_l'):
            pass
        elif (unit == 'per_lm'):
            degree_l = _np.arange(len(mtse))
            mtse /= (2.0 * degree_l + 1.0)
            sd /= (2.0 * degree_l + 1.0)
        else:
            raise ValueError(
                "unit must be 'per_l' or 'per_lm'." +
                "Input value was {:s}".format(repr(unit)))

        if (convention == 'power'):
            return mtse, sd
        elif (convention == 'energy'):
            return mtse * 4.0 * _np.pi, sd * 4.0 * _np.pi
        else:
            raise ValueError(
                "convention must be 'power' or 'energy'." +
                "Input value was {:s}".format(repr(convention)))
1009971
def generate_all_kmers(k):
    
    bases = ['A', 'C', 'G', 'T']
    return [''.join(p) for p in itertools.product(bases, repeat=k)]
262009
def infer_spectra(self, ds):
        
        lvec_all = _get_lvec(ds.test_label_vals, self.pivots, self.scales, derivs=False)
        self.model_spectra = np.dot(lvec_all, self.coeffs.T)
1241492
def eeg_microstates_clustering(data, n_microstates=4, clustering_method="kmeans", n_jobs=1, n_init=25, occurence_rejection_treshold=0.05, max_refitting=5, verbose=True):
    
    # Create training set
    training_set = data.copy()

    if verbose is True:
        print("- Initializing the clustering algorithm...")
    if clustering_method == "kmeans":
        algorithm = sklearn.cluster.KMeans(init='k-means++', n_clusters=n_microstates, n_init=n_init, n_jobs=n_jobs)
    elif clustering_method == "spectral":
        algorithm = sklearn.cluster.SpectralClustering(n_clusters=n_microstates, n_init=n_init, n_jobs=n_jobs)
    elif clustering_method == "agglom":
        algorithm = sklearn.cluster.AgglomerativeClustering(n_clusters=n_microstates, linkage="complete")
    elif clustering_method == "dbscan":
        algorithm = sklearn.cluster.DBSCAN(min_samples=100)
    elif clustering_method == "affinity":
        algorithm = sklearn.cluster.AffinityPropagation(damping=0.5)
    else:
        print("NeuroKit Error: eeg_microstates(): clustering_method must be 'kmeans', 'spectral', 'dbscan', 'affinity' or 'agglom'")


    refitting = 0  # Initialize the number of refittings
    good_fit_achieved = False
    while good_fit_achieved is False:
        good_fit_achieved = True
        if verbose is True:
            print("- Fitting the classifier...")
        # Fit the algorithm
        algorithm.fit(training_set)

        if verbose is True:
            print("- Clustering back the initial data...")
        # Predict the more likely cluster for each observation
        predicted = algorithm.fit_predict(training_set)

        if verbose is True:
            print("- Check for abnormalities...")
        # Check for abnormalities and prune the training set until none found
        occurences = dict(collections.Counter(predicted))
        masks = [np.array([True]*len(training_set))]
        for microstate in occurences:
            # is the frequency of one microstate inferior to a treshold
            if occurences[microstate] < len(data)*occurence_rejection_treshold:
                good_fit_achieved = False
                refitting += 1  # Increment the refitting
                print("NeuroKit Warning: eeg_microstates(): detected some outliers: refitting the classifier (n=" + str(refitting) + ").")
                masks.append(predicted!=microstate)
        mask = np.all(masks, axis=0)
        training_set = training_set[mask]

    return(algorithm)
1348512
protected int interpretLs(LineParser lp, MessageMgr mm){
		String directory = lp.getArgs();
		IOFileFilter fileFilter = new WildcardFileFilter(new String[]{
					"*.ssc"
		});
		DirectoryLoader dl = new CommonsDirectoryWalker(directory, DirectoryFileFilter.INSTANCE, fileFilter);
		if(dl.getLoadErrors().hasErrors()){
			mm.report(dl.getLoadErrors());
			return 1;
		}
		FileSourceList fsl = dl.load();
		if(dl.getLoadErrors().hasErrors()){
			mm.report(dl.getLoadErrors());
			return 1;
		}
		for(FileSource fs : fsl.getSource()){
			//TODO need to adapt to new source return
			mm.report(MessageMgr.createInfoMessage("script file - dir <{}> file <{}>", directory, fs.getBaseFileName()));
		}
		return 0;
	}
1120549
def run_rk4(path,name,E0,laser_frequencies,  N_iter,dt,N_states,
				spectrum_of_laser=None,N_delta=None,frequency_step=None,frequency_end=None,
				rho0=None,print_steps=False,integrate=False,
				save_systems=False):
	""""""

	t0=time()
	params =str(N_iter)+'\n'
	params+=str(dt)+'\n'
	#We give the flag on wether to print each time step.
	if print_steps:
		params+='.true.\n'
	else:
		params+='.false.\n'
	
	#We give the initial value of rho
	N_vars=N_states*(N_states+1)/2-1
	if rho0==None:
		params+=''.join(['(0.0,0.0) ' for i in range(N_vars)])
	elif len(rho0)==N_states-1:
		if sage_included:
			params+=''.join(['('+str(real(i))+','+str(imag(i))+') ' for i in rho0])
		else:
			params+=''.join(['('+str(i.real)+','+str(i.imag)+') ' for i in rho0])
		params+=''.join(['(0.0,0.0) ' for i in range(N_vars-N_states+1)])
	elif len(rho0)==N_vars:
		params+=''.join(['('+str(real(i))+','+str(imag(i))+') ' for i in rho0])

	params+='\n'
	#We give the amplitude of the electrical fields.
	params+=''.join([str(i)+' ' for i in E0])+'\n'
	
	#We give the detuning of each laser (taken from the lowest frequency transition).
	params+=''.join([str(i)+' ' for i in laser_frequencies])+'\n'

	#We give the flag on wether to calculate spectrums or time evolution.
	if spectrum_of_laser==None:
		params+='.false.'
	else:
		if frequency_end !=None:
			if frequency_step !=None:
				raise ValueError,'both frequency_end and frequency_step were specified.'
			if N_delta==1:
				frequency_step=0.0
			else:
				frequency_step=(frequency_end-laser_frequencies[spectrum_of_laser-1])/(N_delta-1)
			
		#frequency_step=frequency_end
		params+='.true.\n'
		params+=str(spectrum_of_laser)+'\n'
		params+=str(N_delta)+'\n'
		params+=str(frequency_step)
	#print params
	
	f=file(path+name+'_params.dat','w')
	f.write(params)
	f.close()
	
	os.system(path+name)
	
	return time()-t0
392608
def find_spectrum_match(spec, spec_lib, method='euclidian'):
    
    # filter out any points with abundance below 1 %
    # spec[spec / np.sum(spec) < 0.01] = 0
    # normalize everything to sum to 1
    spec = spec / np.max(spec)

    if method == 'dot':
        d1 = (spec_lib * lil_matrix(spec).T).sum(axis=1).A ** 2
        d2 = np.sum(spec ** 2) * spec_lib.multiply(spec_lib).sum(axis=1).A
        dist = d1 / d2
    elif method == 'euclidian':
        # st_spc = spectrum[np.newaxis, :].repeat(spec_lib.shape[0], axis=0)
        st_spc = dia_matrix((spec, [0]), shape=(len(spec), len(spec)))
        # calculate the residual sum of squares from spectrum to library
        dist_sp = spec_lib.multiply(spec_lib) - 2 * spec_lib.dot(st_spc)
        dist = dist_sp.sum(axis=1).A + np.sum(spec ** 2)
    return (dist.argmin(), dist.min())
1030025
def do_diagnostics(self):
        ''''''
        self.OLR = self.subprocess['LW'].flux_to_space
        self.LW_down_sfc = self.subprocess['LW'].flux_to_sfc
        self.LW_up_sfc = self.subprocess['LW'].flux_from_sfc
        self.LW_absorbed_sfc = self.LW_down_sfc - self.LW_up_sfc
        self.LW_absorbed_atm = self.subprocess['LW'].absorbed
        self.LW_emission = self.subprocess['LW'].emission
            #  contributions to OLR from surface and atm. levels
            #self.diagnostics['OLR_sfc'] = self.flux['sfc2space']
            #self.diagnostics['OLR_atm'] = self.flux['atm2space']
        self.ASR = (self.subprocess['SW'].flux_from_space -
                    self.subprocess['SW'].flux_to_space)
        #self.SW_absorbed_sfc = (self.subprocess['surface'].SW_from_atm -
        #                            self.subprocess['surface'].SW_to_atm)
        self.SW_absorbed_atm = self.subprocess['SW'].absorbed
        self.SW_down_sfc = self.subprocess['SW'].flux_to_sfc
        self.SW_up_sfc = self.subprocess['SW'].flux_from_sfc
        self.SW_absorbed_sfc = self.SW_down_sfc - self.SW_up_sfc
        self.SW_up_TOA = self.subprocess['SW'].flux_to_space
        self.SW_down_TOA = self.subprocess['SW'].flux_from_space
        self.planetary_albedo = (self.subprocess['SW'].flux_to_space /
                                 self.subprocess['SW'].flux_from_space)
1271916
def average_spectrogram(timeseries, method_func, stride, *args, **kwargs):
    
    # unpack CSD TimeSeries pair, or single timeseries
    try:
        timeseries, other = timeseries
    except ValueError:
        timeseries = timeseries
        other = None

    from ...spectrogram import Spectrogram

    nproc = kwargs.pop('nproc', 1)

    # get params
    epoch = timeseries.t0.value
    nstride = seconds_to_samples(stride, timeseries.sample_rate)
    kwargs['fftlength'] = kwargs.pop('fftlength', stride) or stride
    normalize_fft_params(timeseries, kwargs=kwargs, func=method_func)
    nfft = kwargs['nfft']
    noverlap = kwargs['noverlap']

    # sanity check parameters
    if nstride > timeseries.size:
        raise ValueError("stride cannot be greater than the duration of "
                         "this TimeSeries")
    if nfft > nstride:
        raise ValueError("fftlength cannot be greater than stride")
    if noverlap >= nfft:
        raise ValueError("overlap must be less than fftlength")

    # set up single process Spectrogram method
    def _psd(series):
        """Calculate a single PSD for a spectrogram
        """
        psd_ = _psdn(series, method_func, *args, **kwargs)
        del psd_.epoch  # fixes Segmentation fault (no idea why it faults)
        return psd_

    # define chunks
    tschunks = _chunk_timeseries(timeseries, nstride, noverlap)
    if other is not None:
        otherchunks = _chunk_timeseries(other, nstride, noverlap)
        tschunks = zip(tschunks, otherchunks)

    # calculate PSDs
    psds = mp_utils.multiprocess_with_queues(nproc, _psd, tschunks)

    # recombobulate PSDs into a spectrogram
    return Spectrogram.from_spectra(*psds, epoch=epoch, dt=stride)
1363918
def _multitaper_spectrum(self, clm, k, convention='power', unit='per_l',
                             clat=None, clon=None, coord_degrees=True,
                             lmax=None, taper_wt=None):
        
        if lmax is None:
            lmax = clm.lmax

        if (clat is not None and clon is not None and clat == self.clat and
                clon == self.clon and coord_degrees is self.coord_degrees and
                k <= self.nwinrot):
            # use the already stored coeffs
            pass
        elif (clat is None and clon is None) and \
                (self.clat is not None and self.clon is not None and
                 k <= self.nwinrot):
            # use the already stored coeffs
            pass
        else:
            if clat is None:
                clat = self.clat
            if clon is None:
                clon = self.clon
            if (clat is None and clon is not None) or \
                    (clat is not None and clon is None):
                raise ValueError('clat and clon must both be input. ' +
                                 'clat = {:s}, clon = {:s}'
                                 .format(repr(clat), repr(clon)))
            if clat is None and clon is None:
                self.rotate(clat=90., clon=0., coord_degrees=True, nwinrot=k)
            else:
                self.rotate(clat=clat, clon=clon, coord_degrees=coord_degrees,
                            nwinrot=k)

        sh = clm.to_array(normalization='4pi', csphase=1, lmax=lmax)

        if taper_wt is None:
            mtse, sd = _shtools.SHMultiTaperMaskSE(sh, self.coeffs,
                                                   lmax=lmax, k=k)
        else:
            mtse, sd = _shtools.SHMultiTaperMaskSE(sh, self.coeffs, lmax=lmax,
                                                   k=k, taper_wt=taper_wt)

        if (unit == 'per_l'):
            pass
        elif (unit == 'per_lm'):
            degree_l = _np.arange(len(mtse))
            mtse /= (2.0 * degree_l + 1.0)
            sd /= (2.0 * degree_l + 1.0)
        else:
            raise ValueError(
                "unit must be 'per_l' or 'per_lm'." +
                "Input value was {:s}".format(repr(unit)))

        if (convention == 'power'):
            return mtse, sd
        elif (convention == 'energy'):
            return mtse * 4.0 * _np.pi, sd * 4.0 * _np.pi
        else:
            raise ValueError(
                "convention must be 'power' or 'energy'." +
                "Input value was {:s}".format(repr(convention)))
1254165
def iter_detector_clss():
    
    return iter_subclasses(
        os.path.dirname(os.path.abspath(__file__)),
        Detector,
        _is_abstract_detector,
    )
497785
def main():
    """"""
    parser = argparse.ArgumentParser(description='DistanceClassifier for classification based on distance measure in feature space.',
                                     add_help=False)

    parser.add_argument('INPUT_FILE', type=str, help='Data file to perform DistanceClassifier on; ensure that the class label column is labeled as "class".')

    parser.add_argument('-h', '--help', action='help', help='Show this help message and exit.')

    parser.add_argument('-is', action='store', dest='INPUT_SEPARATOR', default='\t',
                        type=str, help='Character used to separate columns in the input file.')

    parser.add_argument('-d', action='store', dest='D', default='mahalanobis',choices = ['mahalanobis','euclidean'],
                        type=str, help='Distance metric to use.')

    parser.add_argument('-v', action='store', dest='VERBOSITY', default=1, choices=[0, 1, 2],
                        type=int, help='How much information DistanceClassifier communicates while it is running: 0 = none, 1 = minimal, 2 = all.')

    parser.add_argument('-s', action='store', dest='RANDOM_STATE', default=0,
                        type=int, help='Random state for train/test split.')

    parser.add_argument('--version', action='version', version='DistanceClassifier {version}'.format(version=__version__),
                        help='Show DistanceClassifier\'s version number and exit.')

    args = parser.parse_args()

    if args.VERBOSITY >= 2:
        print('\nDistanceClassifier settings:')
        for arg in sorted(args.__dict__):
            print('{}\t=\t{}'.format(arg, args.__dict__[arg]))
        print('')

    input_data = pd.read_csv(args.INPUT_FILE, sep=args.INPUT_SEPARATOR)

    if 'Class' in input_data.columns.values:
        input_data.rename(columns={'Label': 'label'}, inplace=True)

    RANDOM_STATE = args.RANDOM_STATE if args.RANDOM_STATE > 0 else None
    #
    # training_indices, testing_indices = train_test_split(input_data.index,
    #                                                      stratify=input_data['label'].values,
    #                                                      train_size=0.75,
    #                                                      test_size=0.25,
    #                                                      random_state=RANDOM_STATE)
    #
    # training_features = input_data.loc[training_indices].drop('label', axis=1).values
    # training_classes = input_data.loc[training_indices, 'label'].values
    #
    # testing_features = input_data.loc[testing_indices].drop('label', axis=1).values
    # testing_classes = input_data.loc[testing_indices, 'label'].values

    # Run and evaluate DistanceClassifier on the training and testing data
    dc = DistanceClassifier(d = args.D)
    # dc.fit(training_features, training_classes)
    dc.fit(input_data.drop('label',axis=1).values, input_data['label'].values)
    print(dc.score(input_data.drop('label',axis=1).values, input_data['label'].values))
510011
private double ssqerr(int k0, int kmax, double[] logk, double[] log_kDist, double m, double t) {
    int k = kmax - k0;
    double result = 0;
    for(int i = 0; i < k; i++) {
      // double h = log_kDist[i] - (m * (logk[i] - logk[0]) + t); ???
      double h = log_kDist[i] - m * logk[i] - t;
      result += h * h;
    }
    return result;
  }