use of gnu.trove.map.hash.TIntIntHashMap in project mixcr by milaboratory.
the class FullSeqAssembler method calculateRawData.
public RawVariantsData calculateRawData(Supplier<OutputPort<VDJCAlignments>> alignments) {
if (!sequenceToVariantId.isEmpty())
throw new IllegalStateException();
for (byte letter = 0; letter < NucleotideSequence.ALPHABET.basicSize(); letter++) {
NucleotideSequence seq = new NucleotideSequence(new byte[] { letter });
sequenceToVariantId.put(seq, letter);
variantIdToSequence.put(letter, seq);
}
sequenceToVariantId.put(NucleotideSequence.EMPTY, NucleotideSequence.ALPHABET.basicSize());
variantIdToSequence.put(NucleotideSequence.ALPHABET.basicSize(), NucleotideSequence.EMPTY);
TIntIntHashMap coverage = new TIntIntHashMap();
TIntObjectHashMap<TIntObjectHashMap<VariantAggregator>> variants = new TIntObjectHashMap<>();
int nAlignments = 0;
for (VDJCAlignments al : CUtils.it(alignments.get())) {
++nAlignments;
for (PointSequence point : toPointSequences(al)) {
int seqIndex = getVariantIndex(point.sequence.getSequence());
coverage.adjustOrPutValue(point.point, 1, 1);
TIntObjectHashMap<VariantAggregator> map = variants.get(point.point);
if (map == null)
variants.put(point.point, map = new TIntObjectHashMap<>());
VariantAggregator var = map.get(seqIndex);
if (var == null)
map.put(point.point, var = new VariantAggregator());
var.count += 1;
var.sumQuality += 0x7F & point.quality;
}
}
assert nAlignments > 0;
long[] forSort = new long[coverage.size()];
TIntIntIterator iterator = coverage.iterator();
int i = 0;
while (iterator.hasNext()) {
iterator.advance();
forSort[i++] = -((((long) iterator.value()) << 32) | iterator.key());
}
Arrays.sort(forSort);
int[] pointsArray = Arrays.stream(forSort).mapToInt(l -> (int) (-l)).toArray();
TIntIntHashMap revIndex = new TIntIntHashMap();
for (int j = 0; j < pointsArray.length; j++) revIndex.put(pointsArray[j], j);
int[] coverageArray = Arrays.stream(forSort).mapToInt(l -> (int) ((-l) >> 32)).toArray();
int[][] packedData = new int[pointsArray.length][nAlignments];
for (int[] aPackedData : packedData) Arrays.fill(aPackedData, -1);
i = 0;
for (VDJCAlignments al : CUtils.it(alignments.get())) {
for (PointSequence point : toPointSequences(al)) {
int pointIndex = revIndex.get(point.point);
packedData[pointIndex][i] = (sequenceToVariantId.get(point.sequence.getSequence()) << 8) | (0xFF & point.quality);
}
i++;
}
return new RawVariantsData(nAlignments, pointsArray, coverageArray) {
@Override
OutputPort<int[]> createPort() {
return CUtils.asOutputPort(Arrays.asList(packedData));
}
};
}
use of gnu.trove.map.hash.TIntIntHashMap in project mixcr by milaboratory.
the class CloneAssembler method runClustering.
public void runClustering() {
if (clusteredClonesAccumulators != null)
throw new IllegalStateException("Already clustered.");
if (!preClusteringDone)
throw new IllegalStateException("No preclustering done.");
@SuppressWarnings("unchecked") Clustering clustering = new Clustering(cloneList, new SequenceExtractor<CloneAccumulator, NucleotideSequence>() {
@Override
public NucleotideSequence getSequence(CloneAccumulator object) {
return object.getSequence().getConcatenated().getSequence();
}
}, new CloneClusteringStrategy(parameters.getCloneClusteringParameters(), this));
this.progressReporter = clustering;
List<Cluster<CloneAccumulator>> clusters = clustering.performClustering();
clusteredClonesAccumulators = new ArrayList<>(clusters.size());
idMapping = new TIntIntHashMap(cloneList.size());
for (int i = 0; i < clusters.size(); ++i) {
final Cluster<CloneAccumulator> cluster = clusters.get(i);
final CloneAccumulator head = cluster.getHead();
idMapping.put(head.getCloneIndex(), i);
head.setCloneIndex(i);
final int k = ~i;
cluster.processAllChildren(new TObjectProcedure<Cluster<CloneAccumulator>>() {
@Override
public boolean execute(Cluster<CloneAccumulator> object) {
onClustered(head, object.getHead());
if (parameters.isAddReadsCountOnClustering())
head.mergeCounts(object.getHead());
idMapping.put(object.getHead().getCloneIndex(), k);
return true;
}
});
clusteredClonesAccumulators.add(head);
}
this.progressReporter = null;
}
use of gnu.trove.map.hash.TIntIntHashMap in project mixcr by milaboratory.
the class ActionSlice method sliceClnA.
void sliceClnA(ActionHelper helper) throws Exception {
try (ClnAReader reader = new ClnAReader(params.getInputFileName(), VDJCLibraryRegistry.getDefault());
ClnAWriter writer = new ClnAWriter(params.getOutputFileName())) {
// Getting full clone set
CloneSet cloneSet = reader.readCloneSet();
// old clone id -> new clone id
TIntIntHashMap idMapping = new TIntIntHashMap();
long newNumberOfAlignments = 0;
// Creating new cloneset
List<Clone> clones = new ArrayList<>();
int i = 0;
List<OutputPort<VDJCAlignments>> allAlignmentsList = new ArrayList<>();
for (Long cloneId_ : params.ids) {
int cloneId = (int) ((long) cloneId_);
newNumberOfAlignments += reader.numberOfAlignmentsInClone(cloneId);
Clone clone = cloneSet.get(cloneId);
idMapping.put(clone.getId(), i);
clones.add(clone.setId(i).resetParentCloneSet());
OutputPort<VDJCAlignments> als = reader.readAlignmentsOfClone(cloneId);
final int ii = i;
allAlignmentsList.add(() -> {
VDJCAlignments al = als.take();
if (al == null)
return null;
return al.updateCloneIndex(ii);
});
i++;
}
CloneSet newCloneSet = new CloneSet(clones, cloneSet.getUsedGenes(), cloneSet.getAlignedFeatures(), cloneSet.getAlignmentParameters(), cloneSet.getAssemblerParameters());
OutputPort<VDJCAlignments> allAlignmentsPortRaw = new FlatteningOutputPort<>(CUtils.asOutputPort(allAlignmentsList));
AtomicLong idGen = new AtomicLong();
OutputPort<VDJCAlignments> allAlignmentsPort = () -> {
VDJCAlignments al = allAlignmentsPortRaw.take();
if (al == null)
return null;
return al.setAlignmentsIndex(idGen.getAndIncrement());
};
writer.writeClones(newCloneSet);
writer.sortAlignments(allAlignmentsPort, newNumberOfAlignments);
writer.writeAlignmentsAndIndex();
}
}
use of gnu.trove.map.hash.TIntIntHashMap in project BuildCraft by BuildCraft.
the class NbtSquishMapWriter method writeListPacked.
private void writeListPacked(WrittenType type, DataOutput to, NBTTagList list) throws IOException {
profiler.startSection("list_packed");
to.writeByte(NbtSquishConstants.COMPLEX_LIST_PACKED);
profiler.startSection("header");
profiler.startSection("init");
int[] data = new int[list.tagCount()];
TIntIntHashMap indexes = new TIntIntHashMap();
for (int i = 0; i < list.tagCount(); i++) {
profiler.startSection("entry");
profiler.startSection("index");
int index = map.indexOfTag(list.get(i));
profiler.endSection();
data[i] = index;
if (!indexes.increment(index)) {
indexes.put(index, 1);
}
profiler.endSection();
}
// First try to make a simple table
// First sort the indexes into highest count first
profiler.endStartSection("sort");
List<IndexEntry> entries = new ArrayList<>();
for (int index : indexes.keys()) {
int count = indexes.get(index);
IndexEntry entry = new IndexEntry(index, count);
entries.add(entry);
}
entries.sort(Comparator.reverseOrder());
if (debug)
log("\n " + entries.size() + " List entries");
writeVarInt(to, entries.size());
profiler.endStartSection("write");
TIntArrayList sortedIndexes = new TIntArrayList();
int i = 0;
for (IndexEntry entry : entries) {
final int j = i;
NBTBase base = map.getTagForWriting(entry.index);
String n = safeToString(base);
if (debug)
log("\n List entry #" + j + " = " + entry.count + "x" + entry.index + " (" + n + ")");
sortedIndexes.add(entry.index);
type.writeIndex(to, entry.index);
i++;
}
TIntArrayList nextData = new TIntArrayList();
nextData.add(data);
writeVarInt(to, data.length);
profiler.endSection();
profiler.endStartSection("contents");
for (int b = 1; !nextData.isEmpty(); b++) {
profiler.startSection("entry");
CompactingBitSet bitset = new CompactingBitSet(b);
bitset.ensureCapacityValues(nextData.size());
TIntArrayList nextNextData = new TIntArrayList();
int maxVal = (1 << b) - 1;
profiler.startSection("iter");
for (int d : nextData.toArray()) {
// profiler.startSection("entry");
// profiler.startSection("index");
int index = sortedIndexes.indexOf(d);
// profiler.endSection();
if (index < maxVal) {
// profiler.startSection("bitset_append");
bitset.append(index);
// profiler.endSection();
} else {
// profiler.startSection("bitset_append");
bitset.append(maxVal);
// profiler.endStartSection("next_add");
nextNextData.add(d);
// profiler.endSection();
}
// profiler.endSection();
}
profiler.endSection();
sortedIndexes.remove(0, Math.min(sortedIndexes.size(), maxVal));
byte[] bitsetBytes = bitset.getBytes();
if (debug)
log("\n List bitset #" + (bitset.bits - 1));
writeVarInt(to, bitsetBytes.length);
to.write(bitsetBytes);
nextData = nextNextData;
profiler.endSection();
}
profiler.endSection();
profiler.endSection();
}
use of gnu.trove.map.hash.TIntIntHashMap in project GDSC-SMLM by aherbert.
the class CreateData method run.
@Override
public void run(String arg) {
SmlmUsageTracker.recordPlugin(this.getClass(), arg);
extraOptions = ImageJUtils.isExtraOptions();
simpleMode = (arg != null && arg.contains("simple"));
benchmarkMode = (arg != null && arg.contains("benchmark"));
spotMode = (arg != null && arg.contains("spot"));
trackMode = (arg != null && arg.contains("track"));
if ("load".equals(arg)) {
loadBenchmarkData();
return;
}
// Each localisation set is a collection of localisations that represent all localisations
// with the same ID that are on in the same image time frame (Note: the simulation
// can create many localisations per fluorophore per time frame which is useful when
// modelling moving particles)
List<LocalisationModelSet> localisationSets = null;
// Each fluorophore contains the on and off times when light was emitted
List<? extends FluorophoreSequenceModel> fluorophores = null;
if (simpleMode || benchmarkMode || spotMode) {
if (!showSimpleDialog()) {
return;
}
resetMemory();
// 1 second frames
settings.setExposureTime(1000);
areaInUm = settings.getSize() * settings.getPixelPitch() * settings.getSize() * settings.getPixelPitch() / 1e6;
// Number of spots per frame
int count = 0;
int[] nextN = null;
SpatialDistribution dist;
if (benchmarkMode) {
// --------------------
// BENCHMARK SIMULATION
// --------------------
// Draw the same point on the image repeatedly
count = 1;
dist = createFixedDistribution();
try {
reportAndSaveFittingLimits(dist);
} catch (final Exception ex) {
// This will be from the computation of the CRLB
IJ.error(TITLE, ex.getMessage());
return;
}
} else if (spotMode) {
// ---------------
// SPOT SIMULATION
// ---------------
// The spot simulation draws 0 or 1 random point per frame.
// Ensure we have 50% of the frames with a spot.
nextN = new int[settings.getParticles() * 2];
Arrays.fill(nextN, 0, settings.getParticles(), 1);
RandomUtils.shuffle(nextN, UniformRandomProviders.create());
// Only put spots in the central part of the image
final double border = settings.getSize() / 4.0;
dist = createUniformDistribution(border);
} else {
// -----------------
// SIMPLE SIMULATION
// -----------------
// The simple simulation draws n random points per frame to achieve a specified density.
// No points will appear in multiple frames.
// Each point has a random number of photons sampled from a range.
// We can optionally use a mask. Create his first as it updates the areaInUm
dist = createDistribution();
// Randomly sample (i.e. not uniform density in all frames)
if (settings.getSamplePerFrame()) {
final double mean = areaInUm * settings.getDensity();
ImageJUtils.log("Mean samples = %f", mean);
if (mean < 0.5) {
final GenericDialog gd = new GenericDialog(TITLE);
gd.addMessage("The mean samples per frame is low: " + MathUtils.rounded(mean) + "\n \nContinue?");
gd.enableYesNoCancel();
gd.hideCancelButton();
gd.showDialog();
if (!gd.wasOKed()) {
return;
}
}
final PoissonSampler poisson = new PoissonSampler(createRandomGenerator(), mean);
final StoredDataStatistics samples = new StoredDataStatistics(settings.getParticles());
while (samples.getSum() < settings.getParticles()) {
samples.add(poisson.sample());
}
nextN = new int[samples.getN()];
for (int i = 0; i < nextN.length; i++) {
nextN[i] = (int) samples.getValue(i);
}
} else {
// Use the density to get the number per frame
count = (int) Math.max(1, Math.round(areaInUm * settings.getDensity()));
}
}
UniformRandomProvider rng = null;
localisationSets = new ArrayList<>(settings.getParticles());
final int minPhotons = (int) settings.getPhotonsPerSecond();
final int range = (int) settings.getPhotonsPerSecondMaximum() - minPhotons + 1;
if (range > 1) {
rng = createRandomGenerator();
}
// Add frames at the specified density until the number of particles has been reached
int id = 0;
int time = 0;
while (id < settings.getParticles()) {
// Allow the number per frame to be specified
if (nextN != null) {
if (time >= nextN.length) {
break;
}
count = nextN[time];
}
// Simulate random positions in the frame for the specified density
time++;
for (int j = 0; j < count; j++) {
final double[] xyz = dist.next();
// Ignore within border. We do not want to draw things we cannot fit.
// if (!distBorder.isWithinXy(xyz))
// continue;
// Simulate random photons
final int intensity = minPhotons + ((rng != null) ? rng.nextInt(range) : 0);
final LocalisationModel m = new LocalisationModel(id, time, xyz, intensity, LocalisationModel.CONTINUOUS);
// Each localisation can be a separate localisation set
final LocalisationModelSet set = new LocalisationModelSet(id, time);
set.add(m);
localisationSets.add(set);
id++;
}
}
} else {
if (!showDialog()) {
return;
}
resetMemory();
areaInUm = settings.getSize() * settings.getPixelPitch() * settings.getSize() * settings.getPixelPitch() / 1e6;
int totalSteps;
double correlation = 0;
ImageModel imageModel;
if (trackMode) {
// ----------------
// TRACK SIMULATION
// ----------------
// In track mode we create fixed lifetime fluorophores that do not overlap in time.
// This is the simplest simulation to test moving molecules.
settings.setSeconds((int) Math.ceil(settings.getParticles() * (settings.getExposureTime() + settings.getTOn()) / 1000));
totalSteps = 0;
final double simulationStepsPerFrame = (settings.getStepsPerSecond() * settings.getExposureTime()) / 1000.0;
imageModel = new FixedLifetimeImageModel(settings.getStepsPerSecond() * settings.getTOn() / 1000.0, simulationStepsPerFrame, createRandomGenerator());
} else {
// ---------------
// FULL SIMULATION
// ---------------
// The full simulation draws n random points in space.
// The same molecule may appear in multiple frames, move and blink.
//
// Points are modelled as fluorophores that must be activated and then will
// blink and photo-bleach. The molecules may diffuse and this can be simulated
// with many steps per image frame. All steps from a frame are collected
// into a localisation set which can be drawn on the output image.
final SpatialIllumination activationIllumination = createIllumination(settings.getPulseRatio(), settings.getPulseInterval());
// Generate additional frames so that each frame has the set number of simulation steps
totalSteps = (int) Math.ceil(settings.getSeconds() * settings.getStepsPerSecond());
// Since we have an exponential decay of activations
// ensure half of the particles have activated by 30% of the frames.
final double eAct = totalSteps * 0.3 * activationIllumination.getAveragePhotons();
// Q. Does tOn/tOff change depending on the illumination strength?
imageModel = new ActivationEnergyImageModel(eAct, activationIllumination, settings.getStepsPerSecond() * settings.getTOn() / 1000.0, settings.getStepsPerSecond() * settings.getTOffShort() / 1000.0, settings.getStepsPerSecond() * settings.getTOffLong() / 1000.0, settings.getNBlinksShort(), settings.getNBlinksLong(), createRandomGenerator());
imageModel.setUseGeometricDistribution(settings.getNBlinksGeometricDistribution());
// Only use the correlation if selected for the distribution
if (PHOTON_DISTRIBUTION[PHOTON_CORRELATED].equals(settings.getPhotonDistribution())) {
correlation = settings.getCorrelation();
}
}
imageModel.setPhotonBudgetPerFrame(true);
imageModel.setDiffusion2D(settings.getDiffuse2D());
imageModel.setRotation2D(settings.getRotate2D());
IJ.showStatus("Creating molecules ...");
final SpatialDistribution distribution = createDistribution();
final List<CompoundMoleculeModel> compounds = createCompoundMolecules();
if (compounds == null) {
return;
}
final List<CompoundMoleculeModel> molecules = imageModel.createMolecules(compounds, settings.getParticles(), distribution, settings.getRotateInitialOrientation());
// Activate fluorophores
IJ.showStatus("Creating fluorophores ...");
// Note: molecules list will be converted to compounds containing fluorophores
fluorophores = imageModel.createFluorophores(molecules, totalSteps);
if (fluorophores.isEmpty()) {
IJ.error(TITLE, "No fluorophores created");
return;
}
// Map the fluorophore ID to the compound for mixtures
if (compounds.size() > 1) {
idToCompound = new TIntIntHashMap(fluorophores.size());
for (final FluorophoreSequenceModel l : fluorophores) {
idToCompound.put(l.getId(), l.getLabel());
}
}
IJ.showStatus("Creating localisations ...");
// TODO - Output a molecule Id for each fluorophore if using compound molecules. This allows
// analysis
// of the ratio of trimers, dimers, monomers, etc that could be detected.
totalSteps = checkTotalSteps(totalSteps, fluorophores);
if (totalSteps == 0) {
return;
}
imageModel.setPhotonDistribution(createPhotonDistribution());
try {
imageModel.setConfinementDistribution(createConfinementDistribution());
} catch (final ConfigurationException ex) {
// We asked the user if it was OK to continue and they said no
return;
}
// This should be optimised
imageModel.setConfinementAttempts(10);
final List<LocalisationModel> localisations = imageModel.createImage(molecules, settings.getFixedFraction(), totalSteps, settings.getPhotonsPerSecond() / settings.getStepsPerSecond(), correlation, settings.getRotateDuringSimulation());
// Re-adjust the fluorophores to the correct time
if (settings.getStepsPerSecond() != 1) {
final double scale = 1.0 / settings.getStepsPerSecond();
for (final FluorophoreSequenceModel f : fluorophores) {
f.adjustTime(scale);
}
}
// Integrate the frames
localisationSets = combineSimulationSteps(localisations);
localisationSets = filterToImageBounds(localisationSets);
}
datasetNumber.getAndIncrement();
final List<LocalisationModel> localisations = drawImage(localisationSets);
if (localisations == null || localisations.isEmpty()) {
IJ.error(TITLE, "No localisations created");
return;
}
fluorophores = removeFilteredFluorophores(fluorophores, localisations);
final double signalPerFrame = showSummary(fluorophores, localisations);
if (!benchmarkMode) {
final boolean fullSimulation = (!(simpleMode || spotMode));
saveSimulationParameters(localisations.size(), fullSimulation, signalPerFrame);
}
IJ.showStatus("Saving data ...");
saveFluorophores(fluorophores);
saveImageResults(results);
saveLocalisations(localisations);
// The settings for the filenames may have changed
SettingsManager.writeSettings(settings.build());
IJ.showStatus("Done");
}
Aggregations