Computes the edit distance between the two given sequences.
    This uses the relatively fast method that only constructs
    two columns of the 2d array for edits.  This function actually uses four columns
    because we track the number of matches too.
199921
def edit_distance(seq1, seq2, action_function=lowest_cost_action, test=operator.eq):
    
    m = len(seq1)
    n = len(seq2)
    # Special, easy cases:
    if seq1 == seq2:
        return 0, n
    if m == 0:
        return n, 0
    if n == 0:
        return m, 0
    v0 = [0] * (n + 1)     # The two 'error' columns
    v1 = [0] * (n + 1)
    m0 = [0] * (n + 1)     # The two 'match' columns
    m1 = [0] * (n + 1)
    for i in range(1, n + 1):
        v0[i] = i
    for i in range(1, m + 1):
        v1[0] = i
        for j in range(1, n + 1):
            cost = 0 if test(seq1[i - 1], seq2[j - 1]) else 1
            # The costs
            ins_cost = v1[j - 1] + 1
            del_cost = v0[j] + 1
            sub_cost = v0[j - 1] + cost
            # Match counts
            ins_match = m1[j - 1]
            del_match = m0[j]
            sub_match = m0[j - 1] + int(not cost)

            action = action_function(ins_cost, del_cost, sub_cost, ins_match,
                                     del_match, sub_match, cost)

            if action in [EQUAL, REPLACE]:
                v1[j] = sub_cost
                m1[j] = sub_match
            elif action == INSERT:
                v1[j] = ins_cost
                m1[j] = ins_match
            elif action == DELETE:
                v1[j] = del_cost
                m1[j] = del_match
            else:
                raise Exception('Invalid dynamic programming option returned!')
                # Copy the columns over
        for i in range(0, n + 1):
            v0[i] = v1[i]
            m0[i] = m1[i]
    return v1[n], m1[n]
480356
def mutateSequence(seq, distance):
    
    subProb=distance
    inProb=0.05*distance
    deProb=0.05*distance
    contProb=0.9
    l = []
    bases = [ 'A', 'C', 'T', 'G' ]
    i=0
    while i < len(seq):
        if random.random() < subProb:
            l.append(random.choice(bases))
        else:
            l.append(seq[i])
        if random.random() < inProb:
            l += getRandomSequence(_expLength(0, contProb))[1]
        if random.random() < deProb:
            i += int(_expLength(0, contProb))
        i += 1
    return "".join(l)
729965
def _expand(subsequence, sequence, max_l_dist):
    
    # If given a long sub-sequence and relatively small max distance,
    # use a more complex algorithm better optimized for such cases.
    if len(subsequence) > max(max_l_dist * 2, 10):
        return _expand_long(subsequence, sequence, max_l_dist)
    else:
        return _expand_short(subsequence, sequence, max_l_dist)
85435
public static <S extends Sequence<S>> int calculateScore(S seq1, Range seq1Range, Mutations<S> mutations,
                                                             LinearGapAlignmentScoring<S> scoring) {
        if (!mutations.isEmpty() && mutations.getPositionByIndex(0) < seq1Range.getFrom() - 1)
            throw new IllegalArgumentException();

        final AlignmentIteratorForward<S> iterator = new AlignmentIteratorForward<>(mutations, seq1Range);

        int score = 0;

        while (iterator.advance()) {
            final int mut = iterator.getCurrentMutation();
            switch (Mutation.getRawTypeCode(mut)) {
                case RAW_MUTATION_TYPE_SUBSTITUTION:
                    score += scoring.getScore(Mutation.getFrom(mut), Mutation.getTo(mut));
                    break;

                case RAW_MUTATION_TYPE_DELETION:
                case RAW_MUTATION_TYPE_INSERTION:
                    score += scoring.getGapPenalty();
                    break;

                default:
                    byte c = seq1.codeAt(iterator.getSeq1Position());
                    score += scoring.getScore(c, c);
                    break;
            }
        }

        return score;
    }
54473
function calculateDiff() {
    // we need to take 2 pointers
    // and travel progressively between them ...
    var newFirst = [
      [].concat(array1Extract[0]),
      [].concat(array1Extract[1])
    ];
    var newSecond = [
      [].concat(array2Extract[0]),
      [].concat(array2Extract[1])
    ];
    var array1Length = array1Extract[0] ? array1Extract[0].length : 0;
    var array2Length = array2Extract[0] ? array2Extract[0].length : 0;

    var pos1 = 0;
    var pos2 = 0;
    var previous2 = 0;
    while (pos1 < array1Length) {
      var diff = newFirst[0][pos1] - array2Extract[0][pos2];
      if (Math.abs(diff) < widthBottom) { // there is some overlap
        var overlap;
        if (options.trapezoid) {
          overlap = getOverlapTrapezoid(newFirst[0][pos1], newFirst[1][pos1], newSecond[0][pos2], newSecond[1][pos2], widthTop, widthBottom);
        } else {
          overlap = getOverlap(newFirst[0][pos1], newFirst[1][pos1], newSecond[0][pos2], newSecond[1][pos2], widthTop, widthBottom);
        }
        newFirst[1][pos1] -= overlap;
        newSecond[1][pos2] -= overlap;
        if (pos2 < (array2Length - 1)) {
          pos2++;
        } else {
          pos1++;
          pos2 = previous2;
        }
      } else {
        if (diff > 0 && pos2 < (array2Length - 1)) {
          pos2++;
          previous2 = pos2;
        } else {
          pos1++;
          pos2 = previous2;
        }
      }
    }
    return newSecond;
  }
174901
def compare(s1, s2, seg=True, ignore=False, stopwords=False):
    '''
    
    '''
    if s1 == s2: return 1.0
    
    s1_words = []
    s2_words = []

    if seg:
        s1 = [x for x in jieba.cut(s1, cut_all=False, HMM=False)]
        s2 = [x for x in jieba.cut(s2, cut_all=False, HMM=False)]
    else:
        s1 = s1.split()
        s2 = s2.split()

    # check stopwords
    if not stopwords:
        global _stopwords
        for x in s1: 
            if not x in _stopwords:
                s1_words.append(x)
        for x in s2:
            if not x in _stopwords:
                s2_words.append(x)
    else:
        s1_words = s1 
        s2_words = s2

    assert len(s1) > 0 and len(s2) > 0, "The length of s1 and s2 should > 0."
    return _similarity_distance(s1_words, s2_words, ignore)
1124531
int clipLine(Point2D p0, Point2D p1, int lineExtension, double[] segParams,
			double[] boundaryDistances) {
		if (boundaryDistances != null) {
			boundaryDistances[0] = -1.0;
			boundaryDistances[1] = -1.0;
		}

		double[] tOld = new double[2];// LOCALREFCLASS1(ArrayOf(double), int,
										// tOld, 2);
		int modified = 0;

		Point2D delta = new Point2D(p1.x - p0.x, p1.y - p0.y);

		if (delta.x == 0.0 && delta.y == 0.0) // input line degenerates to a
												// point
		{
			segParams[0] = 0.0;
			segParams[1] = 0.0;
			return contains(p0) ? 4 : 0;
		}

		segParams[0] = ((lineExtension & 1) != 0) ? NumberUtils.negativeInf()
				: 0.0;
		segParams[1] = ((lineExtension & 2) != 0) ? NumberUtils.positiveInf()
				: 1.0;
		tOld[0] = segParams[0];
		tOld[1] = segParams[1];

		if (clipLineAuxiliary(delta.x, xmin - p0.x, segParams)
				&& clipLineAuxiliary(-delta.x, p0.x - xmax, segParams)
				&& clipLineAuxiliary(delta.y, ymin - p0.y, segParams)
				&& clipLineAuxiliary(-delta.y, p0.y - ymax, segParams)) {
			if (segParams[1] < tOld[1]) {
				p1.scaleAdd(segParams[1], delta, p0);
				_snapToBoundary(p1); // needed for accuracy
				modified |= 2;

				if (boundaryDistances != null)
					boundaryDistances[1] = _boundaryDistance(p1);
			}
			if (segParams[0] > tOld[0]) {
				p0.scaleAdd(segParams[0], delta, p0);
				_snapToBoundary(p0); // needed for accuracy
				modified |= 1;

				if (boundaryDistances != null)
					boundaryDistances[0] = _boundaryDistance(p0);
			}
		}

		return modified;
	}
554230
def hamming_distance(s1, s2):
    
    if len(s1) != len(s2):
        raise ValueError("Undefined for sequences of unequal length")
    return sum(ch1 != ch2 for ch1, ch2 in zip(s1, s2))
53585
def calculate_columns(sequence):
    
    columns = {}

    for row in sequence:
        for key in row.keys():
            if key not in columns:
                columns[key] = len(key)

            value_length = len(str(row[key]))
            if value_length > columns[key]:
                columns[key] = value_length

    return columns
254458
public float similarity() {
        int match = 0, mismatch = 0;

        AlignmentIteratorForward<S> iterator = forwardIterator();
        while (iterator.advance()) {
            final int mut = iterator.getCurrentMutation();
            if (mut == NON_MUTATION)
                ++match;
            else
                ++mismatch;
        }

        return 1.0f * match / (match + mismatch);
    }
736869
public static float distanceSq( float[] a , float[]b ) {
		float ret = 0;
		for( int i = 0; i < a.length; i++ ) {
			float d = a[i] - b[i];
			ret += d*d;
		}
		return ret;
	}
814452
def calc_2dsplinecoeffs_c(array2d):
    
    #Set up result arrays
    out= copy.copy(array2d)
    out= numpy.require(out,dtype=numpy.float64,requirements=['C','W'])

    #Set up the C code
    ndarrayFlags= ('C_CONTIGUOUS','WRITEABLE')
    interppotential_calc_2dsplinecoeffs= _lib.samples_to_coefficients
    interppotential_calc_2dsplinecoeffs.argtypes= [ndpointer(dtype=numpy.float64,flags=ndarrayFlags),
                                                   ctypes.c_int,
                                                   ctypes.c_int]

    #Run the C code
    interppotential_calc_2dsplinecoeffs(out,out.shape[1],out.shape[0])

    return out
888464
public static int compactDistance(String s1, String s2) {
    if (s1.length() == 0)
      return s2.length();
    if (s2.length() == 0)
      return s1.length();

    // the maximum edit distance there is any point in reporting.
    int maxdist = Math.min(s1.length(), s2.length()) / 2;
    
    // we allocate just one column instead of the entire matrix, in
    // order to save space.  this also enables us to implement the
    // algorithm somewhat faster.  the first cell is always the
    // virtual first row.
    int s1len = s1.length();
    int[] column = new int[s1len + 1];

    // first we need to fill in the initial column. we use a separate
    // loop for this, because in this case our basis for comparison is
    // not the previous column, but a virtual first column.
    int ix2 = 0;
    char ch2 = s2.charAt(ix2);
    column[0] = 1; // virtual first row
    for (int ix1 = 1; ix1 <= s1len; ix1++) {
      int cost = s1.charAt(ix1 - 1) == ch2 ? 0 : 1;

      // Lowest of three: above (column[ix1 - 1]), aboveleft: ix1 - 1,
      // left: ix1. Latter cannot possibly be lowest, so is
      // ignored.
      column[ix1] = Math.min(column[ix1 - 1], ix1 - 1) + cost;
    }

    // okay, now we have an initialized first column, and we can
    // compute the rest of the matrix.
    int above = 0;
    for (ix2 = 1; ix2 < s2.length(); ix2++) {
      ch2 = s2.charAt(ix2);
      above = ix2 + 1; // virtual first row

      int smallest = s1len * 2; // used to implement cutoff
      for (int ix1 = 1; ix1 <= s1len; ix1++) {
        int cost = s1.charAt(ix1 - 1) == ch2 ? 0 : 1;

        // above:     above
        // aboveleft: column[ix1 - 1]
        // left:      column[ix1]
        int value = Math.min(Math.min(above, column[ix1 - 1]), column[ix1]) +
                    cost;
        column[ix1 - 1] = above; // write previous
        above = value;           // keep current
        smallest = Math.min(smallest, value);
      }
      column[s1len] = above;

      // check if we can stop because we'll be going over the max distance
      if (smallest > maxdist)
        return smallest;
    }

    // ok, we're done
    return above;
  }
1171813
void computeColumnSpans() {
        for (int index = 0; index < _generatedMappings().size(); ++index) {
            ParsedMapping mapping = _generatedMappings().get(index);

            // Mappings do not contain a field for the last generated columnt. We
            // can come up with an optimistic estimate, however, by assuming that
            // mappings are contiguous (i.e. given two consecutive mappings, the
            // first mapping ends where the second one starts).
            if (index + 1 < _generatedMappings().size()) {
                ParsedMapping nextMapping = _generatedMappings().get(index + 1);

                if (mapping.generatedLine == nextMapping.generatedLine) {
                    mapping.lastGeneratedColumn = nextMapping.generatedColumn - 1;
                    continue;
                }
            }

            // The last mapping for each line spans the entire line.
            mapping.lastGeneratedColumn = Integer.MAX_VALUE;
        }
    }
921501
def index_of_match_after_additions
      return unless expected_additions > 0

      (1..expected_additions).each do |i|
        next_difference = item_difference(
          old_array[old_index],
          new_array[new_index + i],
          old_index
        )

        return new_index + i if next_difference.empty?
      end

      nil
    end
1025987
def get_tiltplanes(self, sequence):
        '''
        
        '''
        tilting_planes = []
        distance_map = []

        for i in range(1, len(sequence)):
            distance_map.append([ sequence[i], self.virtual_atoms.get_distance( sequence[0], sequence[i] ) ])

        distance_map = sorted(distance_map, key=lambda x: x[1])

        if len(distance_map) == 4:
            # surface edge case
            # semi-octahedron at surface edge has only one tilting plane to consider
            sorted_dist = [i[0] for i in distance_map]
            if distance_map[-1][1] - distance_map[-2][1] < 0.5:
                # 1st case: max diff < 0.5 Angstrom,
                # meaning all distances to reference atom are similar,
                # therefore the reference atom is above the searched plane
                # and the searched plane consists of other atoms
                tilting_planes.append( [ i[0] for i in distance_map ] )
            else:
                # 2nd case: reference atom belongs to the searched plane,
                # procedure needs to be repeated with the next atom as reference atom
                candidates = [sequence[0], sorted_dist[-1]]
                next_distance_map = []
                next_distance_map.append([ sorted_dist[1], self.virtual_atoms.get_distance( sorted_dist[0], sorted_dist[1] ) ])
                next_distance_map.append([ sorted_dist[2], self.virtual_atoms.get_distance( sorted_dist[0], sorted_dist[2] ) ])
                next_distance_map = sorted(next_distance_map, key=lambda x: x[1])
                next_sorted_dist = [i[0] for i in next_distance_map]

                # the next reference atom is taken above the plane (distances are similar)
                if next_distance_map[1][1] - next_distance_map[0][1] < 0.5: candidates.extend([ next_sorted_dist[0], next_sorted_dist[1] ])

                # the next reference atom is taken in the plane (distances are different)
                else: candidates.extend([ sorted_dist[0], next_sorted_dist[1] ])
                tilting_planes.append(candidates)

        elif len(distance_map) == 5:
            # full octahedron case
            # full octahedron has 3 different tilting planes (perpendicular in ideal case)
            sorted_dist = [i[0] for i in distance_map]

            # 1st plane is found as:
            first_plane = sorted_dist[0:4]
            tilting_planes.append(first_plane)
            distance_map_first_plane = []
            for i in range(1, 4):
                distance_map_first_plane.append([ first_plane[i], self.virtual_atoms.get_distance( first_plane[0], first_plane[i] ) ])
            distance_map_first_plane = sorted(distance_map_first_plane, key=lambda x: x[1])
            sorted_first_plane = [i[0] for i in distance_map_first_plane]

            # 2nd and 3rd planes are found as:
            tilting_planes.append([ sequence[0], sorted_dist[4], first_plane[0], sorted_first_plane[2] ])
            tilting_planes.append([ sequence[0], sorted_dist[4], sorted_first_plane[0], sorted_first_plane[1] ])

        # filter planes by Z according to octahedral spatial compound
        filtered = list(filter(lambda x:
            abs(self.virtual_atoms[ x[0] ].z - self.virtual_atoms[ x[1] ].z) < self.OCTAHEDRON_ATOMS_Z_DIFFERENCE and \
            abs(self.virtual_atoms[ x[1] ].z - self.virtual_atoms[ x[2] ].z) < self.OCTAHEDRON_ATOMS_Z_DIFFERENCE and \
            abs(self.virtual_atoms[ x[2] ].z - self.virtual_atoms[ x[3] ].z) < self.OCTAHEDRON_ATOMS_Z_DIFFERENCE,
            tilting_planes
        )) # Py3
        if len(filtered): tilting_planes = filtered

        return tilting_planes
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)