List of usage examples for org.apache.commons.math3.random Well19937c nextInt
public int nextInt(int n) throws IllegalArgumentException
This default implementation is copied from Apache Harmony java.util.Random (r929253).
Implementation notes:
From source file:com.milaboratory.core.alignment.KMapper.java
/** * Performs an alignment for a part of the target sequence. * * <p>This methods is thread-safe and can be concurrently used by several threads if no new sequences added after * its first invocation.</p>//www.j a va 2 s . c om * * @param sequence target sequence * @param from first nucleotide to align (inclusive) * @param to last nucleotide to align (exclusive) * @return a list of hits found in the target sequence */ public KMappingResult align(NucleotideSequence sequence, int from, int to) { ensureBuilt(); ArrayList<KMappingHit> result = new ArrayList<>(); if (to - from <= kValue) return new KMappingResult(null, result); IntArrayList seedPositions = new IntArrayList((to - from) / minDistance + 2); int seedPosition = from; seedPositions.add(seedPosition); Well19937c random = RandomUtil.getThreadLocalRandom(); while ((seedPosition += random.nextInt(maxDistance + 1 - minDistance) + minDistance) < to - kValue) seedPositions.add(seedPosition); //if (seedPositions.get(seedPositions.size() - 1) != to - kValue) seedPositions.add(to - kValue); int[] seeds = new int[seedPositions.size()]; int kmer; IntArrayList[] candidates = new IntArrayList[sequencesInBase]; //Building candidates arrays (seed) int id, offset; for (int i = 0; i < seeds.length; ++i) { kmer = 0; for (int j = seedPositions.get(i); j < seedPositions.get(i) + kValue; ++j) kmer = kmer << 2 | sequence.codeAt(j); seeds[i] = kmer; if (base[kmer].length == 0) continue; for (int record : base[kmer]) { id = record >>> bitsForOffset; offset = record & offsetMask; if (candidates[id] == null) candidates[id] = new IntArrayList(); candidates[id].add((seedPositions.get(i) - offset) << (32 - bitsForOffset) | i); } } //Selecting best candidates (vote) //int resultId = 0; Info info = new Info(); int cFrom, cTo, siFrom, siTo; int maxRawScore = 0, j, i; double preScore; double maxScore = Float.MIN_VALUE; for (i = candidates.length - 1; i >= 0; --i) { //TODO reconsider conditions: if (candidates[i] != null && candidates[i].size() >= ((minAlignmentLength - kValue + 1) / maxDistance) && candidates[i].size() * matchScore >= maxScore * relativeMinScore) { //Sorting (important) candidates[i].sort(); //Calculating best score and offset values getScoreFromOffsets(candidates[i], info); //Theoretical range of target and reference sequence intersection cFrom = max(info.offset, from); cTo = min(info.offset + refLength[i], to) - kValue; //Calculating number of seeds in this range siTo = siFrom = -1; for (j = seedPositions.size() - 1; j >= 0; --j) if ((seedPosition = seedPositions.get(j)) <= cTo) { if (siTo == -1) siTo = j + 1; if (seedPosition < cFrom) { siFrom = j + 1; break; } } //If siFrom not set, first seed is inside the range of target and //reference sequence intersection if (siFrom == -1) siFrom = 0; //Calculating score without penalty preScore = matchScore * info.score; //+ max(siTo - siFrom - info.score, 0) * mismatchScore; //Selecting candidates if (preScore >= absoluteMinScore) { result.add(new KMappingHit(info.offset, i, (float) preScore, siFrom, siTo)); if (maxRawScore < info.score) maxRawScore = info.score; if (maxScore < preScore) maxScore = preScore; } } } int c, seedIndex, seedIndexMask = (0xFFFFFFFF >>> (bitsForOffset)), offsetDelta, currentSeedPosition, prev; float score; KMappingHit hit; maxScore = 0.0; for (j = result.size() - 1; j >= 0; --j) { hit = result.get(j); //Fulfilling the seed positions array //getting seed positions in intersection of target and reference sequences hit.seedOffsets = new int[hit.to - hit.from]; Arrays.fill(hit.seedOffsets, SEED_NOT_FOUND_OFFSET); for (int k = candidates[hit.id].size() - 1; k >= 0; --k) { // offset value | seed index c = candidates[hit.id].get(k); seedIndex = c & seedIndexMask; //filling seed position array with actual positions of seeds inside intersection if (seedIndex >= result.get(j).from && seedIndex < result.get(j).to) { seedIndex -= hit.from; offsetDelta = abs((c >> (32 - bitsForOffset)) - hit.offset); //check if offset difference is less than max allowed and better seed position is found if (offsetDelta <= maxIndels && (hit.seedOffsets[seedIndex] == SEED_NOT_FOUND_OFFSET || abs(hit.seedOffsets[seedIndex] - hit.offset) > offsetDelta)) hit.seedOffsets[seedIndex] = (c >> (32 - bitsForOffset)); } } //Previous seed position prev = -1; c = -1; for (i = hit.seedOffsets.length - 1; i >= 0; --i) if (hit.seedOffsets[i] != SEED_NOT_FOUND_OFFSET) { if (c != -1) //check for situation like: seedsOffset = [25,25,25,25, 28 ,25] //we have to remove such offsets because it's most likely wrong mapping //c - most left index, prev - middle index and i - right most index //but we iterate in reverse direction if (hit.seedOffsets[c] != hit.seedOffsets[prev] && hit.seedOffsets[prev] != hit.seedOffsets[i] && ((hit.seedOffsets[c] < hit.seedOffsets[prev]) != (hit.seedOffsets[prev] < hit.seedOffsets[i]))) { hit.seedOffsets[prev] = SEED_NOT_FOUND_OFFSET; prev = -1; } c = prev; prev = i; } //Calculating score score = 0.0f; for (int off : hit.seedOffsets) if (off != SEED_NOT_FOUND_OFFSET) score += matchScore; else score += mismatchPenalty; //Floating bounds reward if (floatingLeftBound) { prev = -1; for (i = 0; i < hit.seedOffsets.length; ++i) if (hit.seedOffsets[i] != SEED_NOT_FOUND_OFFSET) { if (prev == -1) { prev = i; continue; } //Calculating score gain for deleting first kMer if (matchScore + abs(hit.seedOffsets[i] - hit.seedOffsets[prev]) * offsetShiftPenalty + (i - prev - 1) * mismatchPenalty <= 0.0f) { //Bad kMer hit.seedOffsets[prev] = SEED_NOT_FOUND_OFFSET; prev = i; continue; } score -= prev * mismatchPenalty; break; } } //Floating bounds reward if (floatingRightBound) { prev = -1; for (i = hit.seedOffsets.length - 1; i >= 0; --i) if (hit.seedOffsets[i] != SEED_NOT_FOUND_OFFSET) { if (prev == -1) { prev = i; continue; } //Calculating score gain for deleting first kMer if (matchScore + abs(hit.seedOffsets[i] - hit.seedOffsets[prev]) * offsetShiftPenalty + (prev - 1 - i) * mismatchPenalty <= 0.0f) { //Bad kMer hit.seedOffsets[prev] = SEED_NOT_FOUND_OFFSET; prev = i; continue; } score -= (hit.seedOffsets.length - 1 - prev) * mismatchPenalty; } } c = -1; prev = -1; //totalIndels = 0; for (i = hit.seedOffsets.length - 1; i >= 0; --i) { if (hit.seedOffsets[i] != SEED_NOT_FOUND_OFFSET) { currentSeedPosition = seedPositions.get(i + hit.from) - hit.seedOffsets[i]; if (c == -1) { c = currentSeedPosition; prev = i; } else if (c <= currentSeedPosition) //Removing paradoxical kMer offsets hit.seedOffsets[i] = SEED_NOT_FOUND_OFFSET; else { //totalIndels += abs(hit.seedOffsets[i] - hit.seedOffsets[prev]); score += abs(hit.seedOffsets[i] - hit.seedOffsets[prev]) * offsetShiftPenalty; c = currentSeedPosition; prev = i; } } } hit.score = score; if (score < absoluteMinScore) result.remove(j); if (maxScore < score) maxScore = score; //if (totalIndels > maxIndels * 2) { // result.remove(j); //} } //Removing candidates with score < maxScore * hitsRange maxScore *= relativeMinScore; for (j = result.size() - 1; j >= 0; --j) { hit = result.get(j); if (hit.score < maxScore) result.remove(j); } //Removing seed conflicts return new KMappingResult(seedPositions, result); }