Longest Common Subsequence accuracy measure: calculate accuracy of each prediction as LCS/length
493630
def accuracy_lcs(self, label, pred):
        
        hit = 0.
        total = 0.
        batch_size = label.shape[0]
        for i in range(batch_size):
            l = self._remove_blank(label[i])
            p = []
            for k in range(self.seq_len):
                p.append(np.argmax(pred[k * batch_size + i]))
            p = self.ctc_label(p)
            hit += self._lcs(p, l) * 1.0 / len(l)
            total += 1.0
        assert total == batch_size
        return hit / total
970732
function longestPositiveIncreasingSubsequence(ns, newStart) {
  let seq = [],
    is = [],
    l = -1,
    pre = new Array(ns.length);

  for (let i = newStart, len = ns.length; i < len; i++) {
    let n = ns[i];
    if (n < 0) continue;
    let j = findGreatestIndexLEQ(seq, n);
    if (j !== -1) pre[i] = is[j];
    if (j === l) {
      l++;
      seq[l] = n;
      is[l]  = i;
    } else if (n < seq[j + 1]) {
      seq[j + 1] = n;
      is[j + 1] = i;
    }
  }

  for (let i = is[l]; l >= 0; i = pre[i], l--) {
    seq[l] = i;
  }

  return seq;
}
565618
def seqres_lines(self):
    """"""
  
    lines = []
  
    for chain in self.keys():
      seq = self[chain]
      serNum = 1
      startidx = 0
      while startidx < len(seq):
        endidx = min(startidx+13, len(seq))
        lines += ["SEQRES  %2i %s %4i  %s\n" % (serNum, chain, len(seq), " ".join(seq[startidx:endidx]))]
        serNum += 1
        startidx += 13
    
    return lines
1371184
def _guess_sequence_type_from_string(self, seq):
        ''''''
        # Define expected residues for each sequence type
        aa_chars = ['P','V','L','I','M','F','Y','W','H','K','R','Q','N','E','D','S','X','*']
        aas = set(itertools.chain(aa_chars, [lower(a) for a in aa_chars]))
        na_chars = ['A', 'T', 'G', 'C', 'N', 'U']
        nas = set(itertools.chain(na_chars, [lower(a) for a in na_chars]))

        num_nucleotide = 0
        num_protein = 0
        count = 0
        for residue in seq:
            if residue in nas:
                num_nucleotide += 1
            elif residue in aas:
                num_protein += 1
            else:
                raise Exception('Encountered unexpected character when attempting to guess sequence type: %s' % (residue))
            count += 1
            if count >300: break
        if float(num_protein) / (num_protein+num_nucleotide) > 0.1:
            return self.PROTEIN_SEQUENCE_TYPE
        else:
            return self.NUCLEOTIDE_SEQUENCE_TYPE
368579
def number_aware_alphabetical_cmp(str1, str2):
    

    def flatten_tokens(tokens):
        l = []
        for token in tokens:
            if isinstance(token, str):
                for char in token:
                    l.append(char)
            else:
                assert isinstance(token, float)
                l.append(token)
        return l

    seq1 = flatten_tokens(tokenize_by_number(str1))
    seq2 = flatten_tokens(tokenize_by_number(str2))

    l = min(len(seq1),len(seq2))

    i = 0

    while i < l:
        if seq1[i] < seq2[i]:
            return -1
        elif seq1[i] > seq2[i]:
            return 1
        i += 1

    if len(seq1) < len(seq2):
        return -1
    elif len(seq1) > len(seq2):
        return 1

    return 0
371979
def inferred_var_seqs_plus_flanks(self, ref_seq, flank_length):
        ''''''
        flank_start = max(0, self.POS - flank_length)
        flank_end = min(len(ref_seq) - 1, self.ref_end_pos() + flank_length)
        seqs = [ref_seq[flank_start:self.POS] + self.REF + ref_seq[self.ref_end_pos() + 1: flank_end + 1]]

        for alt in self.ALT:
            seqs.append(ref_seq[flank_start:self.POS] + alt + ref_seq[self.ref_end_pos() + 1: flank_end + 1])

        return flank_start, seqs
817380
def estimate_library_complexity(df, algorithm="RNA-seq"):
    
    DEFAULT_CUTOFFS = {"RNA-seq": (0.25, 0.40)}
    cutoffs = DEFAULT_CUTOFFS[algorithm]
    if len(df) < 5:
        return {"unique_starts_per_read": 'nan',
                "complexity": "NA"}
    model = sm.ols(formula="starts ~ reads", data=df)
    fitted = model.fit()
    slope = fitted.params["reads"]
    if slope <= cutoffs[0]:
        complexity = "LOW"
    elif slope <= cutoffs[1]:
        complexity = "MEDIUM"
    else:
        complexity = "HIGH"

    # for now don't return the complexity flag
    return {"Unique Starts Per Read": float(slope)}
1230486
def consistency_check(text, word_pairs, err, msg, offset=0):
    """"""
    errors = []

    msg = " ".join(msg.split())

    for w in word_pairs:
        matches = [
            [m for m in re.finditer(w[0], text)],
            [m for m in re.finditer(w[1], text)]
        ]

        if len(matches[0]) > 0 and len(matches[1]) > 0:

            idx_minority = len(matches[0]) > len(matches[1])

            for m in matches[idx_minority]:
                errors.append((
                    m.start() + offset,
                    m.end() + offset,
                    err,
                    msg.format(w[~idx_minority], m.group(0)),
                    w[~idx_minority]))

    return errors
1077170
def _lcs(p, l):
        
        # Dynamic Programming Finding LCS
        if len(p) == 0:
            return 0
        P = np.array(list(p)).reshape((1, len(p)))
        L = np.array(list(l)).reshape((len(l), 1))
        M = np.ndarray(shape=(len(P), len(L)), dtype=np.int32)
        for i in range(M.shape[0]):
            for j in range(M.shape[1]):
                up = 0 if i == 0 else M[i-1, j]
                left = 0 if j == 0 else M[i, j-1]

                if i == 0 or j == 0:
                    M[i, j] = max(up, left, M[i, j])
                else:
                    M[i, j] = M[i, j] + M[i - 1, j - 1]
        return M.max()
78700
def _reduce_output(self, outputs, seq_lengths):
        
        batch_size = outputs.shape[0]
        reduced = []
        # Necessary to iterate over batch because of different sequence lengths
        for i in range(batch_size):
            if self.lstm_reduction == "mean":
                # Average over all non-padding reduced
                # Use dim=0 because first dimension disappears after indexing
                reduced.append(outputs[i, : seq_lengths[i], :].mean(dim=0))
            elif self.lstm_reduction == "max":
                # Max-pool over all non-padding reduced
                # Use dim=0 because first dimension disappears after indexing
                reduced.append(outputs[i, : seq_lengths[i], :].max(dim=0)[0])
            elif self.lstm_reduction == "last":
                # Take the last output of the sequence (before padding starts)
                # NOTE: maybe better to take first and last?
                reduced.append(outputs[i, seq_lengths[i] - 1, :])
            elif self.lstm_reduction == "attention":
                reduced.append(self._attention(outputs[i, : seq_lengths[i], :]))
            else:
                msg = (
                    f"Did not recognize lstm kwarg 'lstm_reduction' == "
                    f"{self.lstm_reduction}"
                )
                raise ValueError(msg)
        return torch.stack(reduced, dim=0)
73950
def _seqfeature_to_coral(feature):
    '''

    '''
    # Some genomic sequences don't have a label attribute
    # TODO: handle genomic cases differently than others. Some features lack
    # a label but should still be incorporated somehow.
    qualifiers = feature.qualifiers
    if 'label' in qualifiers:
        feature_name = qualifiers['label'][0]
    elif 'locus_tag' in qualifiers:
        feature_name = qualifiers['locus_tag'][0]
    else:
        raise FeatureNameError('Unrecognized feature name')
    # Features with gaps are special, require looking at subfeatures
    # Assumption: subfeatures are never more than one level deep
    if feature.location_operator == 'join':
        # Feature has gaps. Have to figure out start/stop from subfeatures,
        # calculate gap indices. A nested feature model may be required
        # eventually.
        # Reorder the sub_feature list by start location
        # Assumption: none of the subfeatures overlap so the last entry in
        # the reordered list also has the final stop point of the feature.
        # FIXME: Getting a deprecation warning about using sub_features
        # instead of feature.location being a CompoundFeatureLocation
        reordered = sorted(feature.location.parts,
                           key=lambda location: location.start)
        starts = [int(location.start) for location in reordered]
        stops = [int(location.end) for location in reordered]
        feature_start = starts.pop(0)
        feature_stop = stops.pop(-1)
        starts = [start - feature_start for start in starts]
        stops = [stop - feature_start for stop in stops]
        feature_gaps = list(zip(stops, starts))
    else:
        # Feature doesn't have gaps. Ignore subfeatures.
        feature_start = int(feature.location.start)
        feature_stop = int(feature.location.end)
        feature_gaps = []
    feature_type = _process_feature_type(feature.type)
    if feature.location.strand == -1:
        feature_strand = 1
    else:
        feature_strand = 0
    if 'gene' in qualifiers:
        gene = qualifiers['gene']
    else:
        gene = []
    if 'locus_tag' in qualifiers:
        locus_tag = qualifiers['locus_tag']
    else:
        locus_tag = []
    coral_feature = coral.Feature(feature_name, feature_start,
                                  feature_stop, feature_type,
                                  gene=gene, locus_tag=locus_tag,
                                  qualifiers=qualifiers,
                                  strand=feature_strand,
                                  gaps=feature_gaps)
    return coral_feature
981747
def pairwise_compare(afa, leven, threads, print_list, ignore_gaps):
    
    # load sequences into dictionary
    seqs = {seq[0]: seq for seq in nr_fasta([afa], append_index = True)}
    num_seqs = len(seqs)
    # define all pairs
    pairs = ((i[0], i[1], ignore_gaps) for i in itertools.combinations(list(seqs.values()), 2))
    pool = multithread(threads)
    # calc percent identity between all pairs - parallelize
    if leven is True:
        pident = pool.map(compare_seqs_leven, pairs)
    else:
        compare = pool.imap_unordered(compare_seqs, pairs)
        pident = [i for i in tqdm(compare, total = (num_seqs*num_seqs)/2)]
    pool.close()
    pool.terminate()
    pool.join()
    return to_dictionary(pident, print_list)
1183189
public static int countMatches(CharSequence sourceSequence, CharSequence subSequence) {
        if (isEmpty(sourceSequence) || isEmpty(subSequence) || sourceSequence.length() < subSequence.length()) {
            return 0;
        }

        int count = 0;
        int sourceSequenceIndex = 0;
        int subSequenceIndex = 0;

        while(sourceSequenceIndex < sourceSequence.length()) {
            if(sourceSequence.charAt(sourceSequenceIndex) == subSequence.charAt(subSequenceIndex)) {
                sourceSequenceIndex++;
                subSequenceIndex++;
                while(sourceSequenceIndex < sourceSequence.length() && subSequenceIndex < subSequence.length()) {
                    if(sourceSequence.charAt(sourceSequenceIndex) != subSequence.charAt(subSequenceIndex)) {
                        break;
                    }
                    sourceSequenceIndex++;
                    subSequenceIndex++;
                }
                if(subSequenceIndex == subSequence.length()) {
                    count++;
                }
                subSequenceIndex = 0;
                continue;
            }
            sourceSequenceIndex++;
        }

        return count;
    }
821658
def _calculate_similarity(c):
    
    ma = {}
    for idc in c:
        set1 = _get_seqs(c[idc])
        [ma.update({(idc, idc2): _common(set1, _get_seqs(c[idc2]), idc, idc2)}) for idc2 in c if idc != idc2 and (idc2, idc) not in ma]
    # logger.debug("_calculate_similarity_ %s" % ma)
    return ma
924318
def max_sequence_length(self, dataset_split):
    
    return {
        problem.DatasetSplit.TRAIN: 64,
        problem.DatasetSplit.EVAL: 128,
        problem.DatasetSplit.TEST: 128
    }[dataset_split]
747399
def plotAccuracyDuringSequenceInference(dirName, title="", yaxis=""):
  
  # Read in results file
  with open(os.path.join(dirName,
            "sequence_batch_high_dec_normal_features.pkl"), "rb") as f:
    results = cPickle.load(f)

  locationRange = []
  featureRange = []
  for r in results:
    if r["numLocations"] not in locationRange: locationRange.append(r["numLocations"])
    if r["numFeatures"] not in featureRange: featureRange.append(r["numFeatures"])
    locationRange.sort()
  featureRange.sort()
  if 10 in featureRange: featureRange.remove(10)
  print "locationRange=",locationRange
  print "featureRange=",featureRange

  ########################################################################
  #
  # Accumulate the L2 accuracies for each condition in a list and compute mean
  # and stdeviations
  # For TM we average across all feature ranges
  L2Accuracies = defaultdict(list)
  TMAccuracies = defaultdict(list)
  for r in results:
    if r["numFeatures"] in featureRange:
      L2Accuracies[(r["numLocations"], r["numFeatures"])].append(r["sequenceAccuracyL2"])
      TMAccuracies[r["numLocations"]].append(r["sequenceCorrectSparsityTM"])

  # meanAccuracy[o,f] = accuracy of TM with o objects and f unique features.
  meanL2Accuracy = numpy.zeros((max(locationRange)+1, max(featureRange) + 1))
  stdevL2 = numpy.zeros((max(locationRange)+1, max(featureRange) + 1))
  meanTMAccuracy = numpy.zeros(max(locationRange)+1)
  stdevTM = numpy.zeros(max(locationRange)+1)
  for o in locationRange:
    for f in featureRange:
      a = numpy.array(L2Accuracies[(o, f)])
      meanL2Accuracy[o, f] = 100.0*a.mean()
      stdevL2[o, f] = 100.0*a.std()

    # Accuracies for TM
    a = numpy.array(TMAccuracies[o])
    meanTMAccuracy[o] = 100.0*a.mean()
    stdevTM[o] = 100.0*a.std()

  ########################################################################
  #
  # Create the plot.
  plt.figure()
  plotPath = os.path.join(os.path.dirname(os.path.realpath(__file__)),
                          "plots", "accuracy_during_sequence_inference.pdf")

  # Plot each curve
  legendList = []
  colorList = ['r', 'b', 'g', 'm', 'c', 'k', 'y']

  for i in range(len(featureRange)):
    f = featureRange[i]
    legendList.append('Sensorimotor layer, feature pool size: {}'.format(f))
    plt.errorbar(locationRange, meanL2Accuracy[locationRange, f],
             yerr=stdevL2[locationRange, f],
             color=colorList[i])

  plt.errorbar(locationRange, meanTMAccuracy[locationRange],
               yerr=stdevTM[locationRange],
               color=colorList[len(featureRange)])
  legendList.append('Temporal sequence layer')

  # format
  plt.legend(legendList, bbox_to_anchor=(0., 0.65, 1., .102), loc="right", prop={'size':10})
  plt.xlabel("Size of location pool")
  # plt.xticks(range(0,max(locationRange)+1,10))
  # plt.yticks(range(0,int(accuracy.max())+2,10))
  plt.ylim(-10.0, 110.0)
  plt.ylabel(yaxis)
  plt.title(title)

    # save
  plt.savefig(plotPath)
  plt.close()
112048
def compare_seqs_leven(seqs):
    
    A, B, ignore_gaps = seqs
    a, b = remove_gaps(A[1], B[1]) # actual sequences
    if len(a) != len(b):
        print('# reads are not the same length', file=sys.stderr)
        exit()
    pident = lr(a, b) * 100
    return A[0], B[0], pident