src/jalview/datamodel/features/FeatureStore.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.datamodel.features;
  22
  23 import jalview.datamodel.SequenceFeature;
  24
  25 import java.util.ArrayList;
  26 import java.util.Arrays;
  27 import java.util.BitSet;
  28 import java.util.Collections;
  29 import java.util.Comparator;
  30 import java.util.HashSet;
  31 import java.util.List;
  32 import java.util.Set;
  33
  34 import intervalstore.api.IntervalStoreI;
  35 import intervalstore.impl.BinarySearcher;
  36 import intervalstore.impl.IntervalStore;
  37
  38
  39 /**
  40  * A data store for a set of sequence features that supports efficient lookup of
  41  * features overlapping a given range. Intended for (but not limited to) storage
  42  * of features for one sequence and feature type.
  43  *
  44  * @author gmcarstairs
  45  *
  46  */
  47 public class FeatureStore
  48 {
  49   /*
  50    * Non-positional features have no (zero) start/end position.
  51    * Kept as a separate list in case this criterion changes in future.
  52    */
  53   List<SequenceFeature> nonPositionalFeatures;
  54
  55   /*
  56    * contact features ordered by first contact position
  57    */
  58   List<SequenceFeature> contactFeatureStarts;
  59
  60   /*
  61    * contact features ordered by second contact position
  62    */
  63   List<SequenceFeature> contactFeatureEnds;
  64
  65   /*
  66    * IntervalStore holds remaining features and provides efficient
  67    * query for features overlapping any given interval
  68    */
  69   IntervalStoreI<SequenceFeature> features;
  70
  71   /*
  72    * Feature groups represented in stored positional features
  73    * (possibly including null)
  74    */
  75   Set<String> positionalFeatureGroups;
  76
  77   /*
  78    * Feature groups represented in stored non-positional features
  79    * (possibly including null)
  80    */
  81   Set<String> nonPositionalFeatureGroups;
  82
  83   /*
  84    * the total length of all positional features; contact features count 1 to
  85    * the total and 1 to size(), consistent with an average 'feature length' of 1
  86    */
  87   int totalExtent;
  88
  89   float positionalMinScore;
  90
  91   float positionalMaxScore;
  92
  93   float nonPositionalMinScore;
  94
  95   float nonPositionalMaxScore;
  96
  97   private ArrayList<SequenceFeature> featuresList;
  98
  99   /**
 100    * Constructor
 101    */
 102   public FeatureStore()
 103   {
 104     features = new IntervalStore<>();
 105     featuresList = new ArrayList<>();
 106     positionalFeatureGroups = new HashSet<>();
 107     nonPositionalFeatureGroups = new HashSet<>();
 108     positionalMinScore = Float.NaN;
 109     positionalMaxScore = Float.NaN;
 110     nonPositionalMinScore = Float.NaN;
 111     nonPositionalMaxScore = Float.NaN;
 112
 113     // we only construct nonPositionalFeatures, contactFeatures if we need to
 114   }
 115
 116   /**
 117    * Adds one sequence feature to the store, and returns true, unless the
 118    * feature is already contained in the store, in which case this method
 119    * returns false. Containment is determined by SequenceFeature.equals()
 120    * comparison.
 121    *
 122    * @param feature
 123    */
 124
 125   public boolean addFeature(SequenceFeature feature)
 126   {
 127     if (contains(feature))
 128     {
 129       return false;
 130     }
 131
 132     /*
 133      * keep a record of feature groups
 134      */
 135     if (!feature.isNonPositional())
 136     {
 137       positionalFeatureGroups.add(feature.getFeatureGroup());
 138     }
 139
 140     if (feature.isContactFeature())
 141     {
 142       addContactFeature(feature);
 143     }
 144     else if (feature.isNonPositional())
 145     {
 146       addNonPositionalFeature(feature);
 147     }
 148     else
 149     {
 150       addNestedFeature(feature);
 151     }
 152
 153     /*
 154      * record the total extent of positional features, to make
 155      * getTotalFeatureLength possible; we count the length of a
 156      * contact feature as 1
 157      */
 158     totalExtent += getFeatureLength(feature);
 159
 160     /*
 161      * record the minimum and maximum score for positional
 162      * and non-positional features
 163      */
 164     float score = feature.getScore();
 165     if (!Float.isNaN(score))
 166     {
 167       if (feature.isNonPositional())
 168       {
 169         nonPositionalMinScore = min(nonPositionalMinScore, score);
 170         nonPositionalMaxScore = max(nonPositionalMaxScore, score);
 171       }
 172       else
 173       {
 174         positionalMinScore = min(positionalMinScore, score);
 175         positionalMaxScore = max(positionalMaxScore, score);
 176       }
 177     }
 178
 179     return true;
 180   }
 181
 182   /**
 183    * Answers true if this store contains the given feature (testing by
 184    * SequenceFeature.equals), else false
 185    *
 186    * @param feature
 187    * @return
 188    */
 189   public boolean contains(SequenceFeature feature)
 190   {
 191     if (feature.isNonPositional())
 192     {
 193       return nonPositionalFeatures == null ? false
 194               : nonPositionalFeatures.contains(feature);
 195     }
 196
 197     if (feature.isContactFeature())
 198     {
 199       return contactFeatureStarts == null ? false
 200               : listContains(contactFeatureStarts, feature);
 201     }
 202
 203     return features == null ? false : features.contains(feature);
 204   }
 205
 206   /**
 207    * Answers the 'length' of the feature, counting 0 for non-positional features
 208    * and 1 for contact features
 209    *
 210    * @param feature
 211    * @return
 212    */
 213   protected static int getFeatureLength(SequenceFeature feature)
 214   {
 215     if (feature.isNonPositional())
 216     {
 217       return 0;
 218     }
 219     if (feature.isContactFeature())
 220     {
 221       return 1;
 222     }
 223     return 1 + feature.getEnd() - feature.getBegin();
 224   }
 225
 226   /**
 227    * Adds the feature to the list of non-positional features (with lazy
 228    * instantiation of the list if it is null), and returns true. The feature
 229    * group is added to the set of distinct feature groups for non-positional
 230    * features. This method allows duplicate features, so test before calling to
 231    * prevent this.
 232    *
 233    * @param feature
 234    */
 235   protected boolean addNonPositionalFeature(SequenceFeature feature)
 236   {
 237     if (nonPositionalFeatures == null)
 238     {
 239       nonPositionalFeatures = new ArrayList<>();
 240     }
 241
 242     nonPositionalFeatures.add(feature);
 243
 244     nonPositionalFeatureGroups.add(feature.getFeatureGroup());
 245
 246     return true;
 247   }
 248
 249   /**
 250    * Adds one feature to the IntervalStore that can manage nested features
 251    * (creating the IntervalStore if necessary)
 252    */
 253   protected synchronized void addNestedFeature(SequenceFeature feature)
 254   {
 255     if (features == null)
 256     {
 257       features = new IntervalStore<>();
 258     }
 259     features.add(feature);
 260     featuresList.add(feature);
 261   }
 262
 263   /**
 264    * Add a contact feature to the lists that hold them ordered by start (first
 265    * contact) and by end (second contact) position, ensuring the lists remain
 266    * ordered, and returns true. This method allows duplicate features to be
 267    * added, so test before calling to avoid this.
 268    *
 269    * @param feature
 270    * @return
 271    */
 272   protected synchronized boolean addContactFeature(SequenceFeature feature)
 273   {
 274     if (contactFeatureStarts == null)
 275     {
 276       contactFeatureStarts = new ArrayList<>();
 277     }
 278     if (contactFeatureEnds == null)
 279     {
 280       contactFeatureEnds = new ArrayList<>();
 281     }
 282
 283     /*
 284      * insert into list sorted by start (first contact position):
 285      * binary search the sorted list to find the insertion point
 286      */
 287     int insertPosition = BinarySearcher.findFirst(contactFeatureStarts,
 288             f -> f.getBegin() >= feature.getBegin());
 289     contactFeatureStarts.add(insertPosition, feature);
 290
 291     /*
 292      * insert into list sorted by end (second contact position):
 293      * binary search the sorted list to find the insertion point
 294      */
 295     insertPosition = BinarySearcher.findFirst(contactFeatureEnds,
 296             f -> f.getEnd() >= feature.getEnd());
 297     contactFeatureEnds.add(insertPosition, feature);
 298
 299     return true;
 300   }
 301
 302   /**
 303    * Answers true if the list contains the feature, else false. This method is
 304    * optimised for the condition that the list is sorted on feature start
 305    * position ascending, and will give unreliable results if this does not hold.
 306    *
 307    * @param features
 308    * @param feature
 309    * @return
 310    */
 311   protected static boolean listContains(List<SequenceFeature> features,
 312           SequenceFeature feature)
 313   {
 314     if (features == null || feature == null)
 315     {
 316       return false;
 317     }
 318
 319     /*
 320      * locate the first entry in the list which does not precede the feature
 321      */
 322     // int pos = binarySearch(features,
 323     // SearchCriterion.byFeature(feature, RangeComparator.BY_START_POSITION));
 324     int pos = BinarySearcher.findFirst(features,
 325             val -> val.getBegin() >= feature.getBegin());
 326     int len = features.size();
 327     while (pos < len)
 328     {
 329       SequenceFeature sf = features.get(pos);
 330       if (sf.getBegin() > feature.getBegin())
 331       {
 332         return false; // no match found
 333       }
 334       if (sf.equals(feature))
 335       {
 336         return true;
 337       }
 338       pos++;
 339     }
 340     return false;
 341   }
 342
 343   /**
 344    * Returns a (possibly empty) list of features whose extent overlaps the given
 345    * range. The returned list is not ordered. Contact features are included if
 346    * either of the contact points lies within the range.
 347    *
 348    * @param start
 349    *          start position of overlap range (inclusive)
 350    * @param end
 351    *          end position of overlap range (inclusive)
 352    * @return
 353    */
 354
 355   public List<SequenceFeature> findOverlappingFeatures(long start, long end)
 356   {
 357     List<SequenceFeature> result = new ArrayList<>();
 358
 359     findContactFeatures(start, end, result);
 360
 361     if (features != null)
 362     {
 363       result.addAll(features.findOverlaps(start, end));
 364     }
 365
 366     return result;
 367   }
 368
 369   /**
 370    * Adds contact features to the result list where either the second or the
 371    * first contact position lies within the target range
 372    *
 373    * @param from
 374    * @param to
 375    * @param result
 376    */
 377   protected void findContactFeatures(long from, long to,
 378           List<SequenceFeature> result)
 379   {
 380     if (contactFeatureStarts != null)
 381     {
 382       findContactStartOverlaps(from, to, result);
 383     }
 384     if (contactFeatureEnds != null)
 385     {
 386       findContactEndOverlaps(from, to, result);
 387     }
 388   }
 389
 390   /**
 391    * Adds to the result list any contact features whose end (second contact
 392    * point), but not start (first contact point), lies in the query from-to
 393    * range
 394    *
 395    * @param from
 396    * @param to
 397    * @param result
 398    */
 399   protected void findContactEndOverlaps(long from, long to,
 400           List<SequenceFeature> result)
 401   {
 402     /*
 403      * find the first contact feature (if any)
 404      * whose end point is not before the target range
 405      */
 406     int index = BinarySearcher.findFirst(contactFeatureEnds,
 407             f -> f.getEnd() >= from);
 408
 409     while (index < contactFeatureEnds.size())
 410     {
 411       SequenceFeature sf = contactFeatureEnds.get(index);
 412       if (!sf.isContactFeature())
 413       {
 414         System.err.println("Error! non-contact feature type " + sf.getType()
 415                 + " in contact features list");
 416         index++;
 417         continue;
 418       }
 419
 420       int begin = sf.getBegin();
 421       if (begin >= from && begin <= to)
 422       {
 423         /*
 424          * this feature's first contact position lies in the search range
 425          * so we don't include it in results a second time
 426          */
 427         index++;
 428         continue;
 429       }
 430
 431       if (sf.getEnd() > to)
 432       {
 433         /*
 434          * this feature (and all following) has end point after the target range
 435          */
 436         break;
 437       }
 438
 439       /*
 440        * feature has end >= from and end <= to
 441        * i.e. contact end point lies within overlap search range
 442        */
 443       result.add(sf);
 444       index++;
 445     }
 446   }
 447
 448   /**
 449    * Adds contact features whose start position lies in the from-to range to the
 450    * result list
 451    *
 452    * @param from
 453    * @param to
 454    * @param result
 455    */
 456   protected void findContactStartOverlaps(long from, long to,
 457           List<SequenceFeature> result)
 458   {
 459     int index = BinarySearcher.findFirst(contactFeatureStarts,
 460             f -> f.getBegin() >= from);
 461
 462     while (index < contactFeatureStarts.size())
 463     {
 464       SequenceFeature sf = contactFeatureStarts.get(index);
 465       if (!sf.isContactFeature())
 466       {
 467         System.err.println("Error! non-contact feature " + sf.toString()
 468                 + " in contact features list");
 469         index++;
 470         continue;
 471       }
 472       if (sf.getBegin() > to)
 473       {
 474         /*
 475          * this feature's start (and all following) follows the target range
 476          */
 477         break;
 478       }
 479
 480       /*
 481        * feature has begin >= from and begin <= to
 482        * i.e. contact start point lies within overlap search range
 483        */
 484       result.add(sf);
 485       index++;
 486     }
 487   }
 488
 489   /**
 490    * Answers a list of all positional features stored, in no guaranteed order
 491    *
 492    * @return
 493    */
 494
 495   public List<SequenceFeature> getPositionalFeatures()
 496   {
 497     List<SequenceFeature> result = new ArrayList<>();
 498
 499     /*
 500      * add any contact features - from the list by start position
 501      */
 502     if (contactFeatureStarts != null)
 503     {
 504       result.addAll(contactFeatureStarts);
 505     }
 506
 507     /*
 508      * add any nested features
 509      */
 510     if (features != null)
 511     {
 512       result.addAll(features);
 513     }
 514
 515     return result;
 516   }
 517
 518   /**
 519    * Answers a list of all contact features. If there are none, returns an
 520    * immutable empty list.
 521    *
 522    * @return
 523    */
 524
 525   public List<SequenceFeature> getContactFeatures()
 526   {
 527     if (contactFeatureStarts == null)
 528     {
 529       return Collections.emptyList();
 530     }
 531     return new ArrayList<>(contactFeatureStarts);
 532   }
 533
 534   /**
 535    * Answers a list of all non-positional features. If there are none, returns
 536    * an immutable empty list.
 537    *
 538    * @return
 539    */
 540
 541   public List<SequenceFeature> getNonPositionalFeatures()
 542   {
 543     if (nonPositionalFeatures == null)
 544     {
 545       return Collections.emptyList();
 546     }
 547     return new ArrayList<>(nonPositionalFeatures);
 548   }
 549
 550   /**
 551    * Deletes the given feature from the store, returning true if it was found
 552    * (and deleted), else false. This method makes no assumption that the feature
 553    * is in the 'expected' place in the store, in case it has been modified since
 554    * it was added.
 555    *
 556    * @param sf
 557    */
 558
 559   public synchronized boolean delete(SequenceFeature sf)
 560   {
 561     boolean removed = false;
 562
 563     /*
 564      * try contact positions (and if found, delete
 565      * from both lists of contact positions)
 566      */
 567     if (!removed && contactFeatureStarts != null)
 568     {
 569       removed = contactFeatureStarts.remove(sf);
 570       if (removed)
 571       {
 572         contactFeatureEnds.remove(sf);
 573       }
 574     }
 575
 576     boolean removedNonPositional = false;
 577
 578     /*
 579      * if not found, try non-positional features
 580      */
 581     if (!removed && nonPositionalFeatures != null)
 582     {
 583       removedNonPositional = nonPositionalFeatures.remove(sf);
 584       removed = removedNonPositional;
 585     }
 586
 587     /*
 588      * if not found, try nested features
 589      */
 590     if (!removed && features != null)
 591     {
 592       removed = features.remove(sf);
 593       featuresList.remove(sf);
 594     }
 595
 596     if (removed)
 597     {
 598       rescanAfterDelete();
 599     }
 600
 601     return removed;
 602   }
 603
 604   /**
 605    * Rescan all features to recompute any cached values after an entry has been
 606    * deleted. This is expected to be an infrequent event, so performance here is
 607    * not critical.
 608    */
 609   protected synchronized void rescanAfterDelete()
 610   {
 611     positionalFeatureGroups.clear();
 612     nonPositionalFeatureGroups.clear();
 613     totalExtent = 0;
 614     positionalMinScore = Float.NaN;
 615     positionalMaxScore = Float.NaN;
 616     nonPositionalMinScore = Float.NaN;
 617     nonPositionalMaxScore = Float.NaN;
 618     /*
 619      * scan non-positional features for groups and scores
 620      */
 621     for (SequenceFeature sf : getNonPositionalFeatures())
 622     {
 623       nonPositionalFeatureGroups.add(sf.getFeatureGroup());
 624       float score = sf.getScore();
 625       nonPositionalMinScore = min(nonPositionalMinScore, score);
 626       nonPositionalMaxScore = max(nonPositionalMaxScore, score);
 627     }
 628
 629     /*
 630      * scan positional features for groups, scores and extents
 631      */
 632     for (SequenceFeature sf : getPositionalFeatures())
 633     {
 634       positionalFeatureGroups.add(sf.getFeatureGroup());
 635       float score = sf.getScore();
 636       positionalMinScore = min(positionalMinScore, score);
 637       positionalMaxScore = max(positionalMaxScore, score);
 638       totalExtent += getFeatureLength(sf);
 639     }
 640   }
 641
 642   /**
 643    * A helper method to return the minimum of two floats, where a non-NaN value
 644    * is treated as 'less than' a NaN value (unlike Math.min which does the
 645    * opposite)
 646    *
 647    * @param f1
 648    * @param f2
 649    */
 650   protected static float min(float f1, float f2)
 651   {
 652     if (Float.isNaN(f1))
 653     {
 654       return Float.isNaN(f2) ? f1 : f2;
 655     }
 656     else
 657     {
 658       return Float.isNaN(f2) ? f1 : Math.min(f1, f2);
 659     }
 660   }
 661
 662   /**
 663    * A helper method to return the maximum of two floats, where a non-NaN value
 664    * is treated as 'greater than' a NaN value (unlike Math.max which does the
 665    * opposite)
 666    *
 667    * @param f1
 668    * @param f2
 669    */
 670   protected static float max(float f1, float f2)
 671   {
 672     if (Float.isNaN(f1))
 673     {
 674       return Float.isNaN(f2) ? f1 : f2;
 675     }
 676     else
 677     {
 678       return Float.isNaN(f2) ? f1 : Math.max(f1, f2);
 679     }
 680   }
 681
 682   /**
 683    * Answers true if this store has no features, else false
 684    *
 685    * @return
 686    */
 687
 688   public boolean isEmpty()
 689   {
 690     boolean hasFeatures = (contactFeatureStarts != null
 691             && !contactFeatureStarts.isEmpty())
 692             || (nonPositionalFeatures != null
 693                     && !nonPositionalFeatures.isEmpty())
 694             || (features != null && features.size() > 0);
 695
 696     return !hasFeatures;
 697   }
 698
 699   /**
 700    * Answers the set of distinct feature groups stored, possibly including null,
 701    * as an unmodifiable view of the set. The parameter determines whether the
 702    * groups for positional or for non-positional features are returned.
 703    *
 704    * @param positionalFeatures
 705    * @return
 706    */
 707
 708   public Set<String> getFeatureGroups(boolean positionalFeatures)
 709   {
 710     if (positionalFeatures)
 711     {
 712       return Collections.unmodifiableSet(positionalFeatureGroups);
 713     }
 714     else
 715     {
 716       return nonPositionalFeatureGroups == null
 717               ? Collections.<String> emptySet()
 718               : Collections.unmodifiableSet(nonPositionalFeatureGroups);
 719     }
 720   }
 721
 722   /**
 723    * Answers the number of positional (or non-positional) features stored.
 724    * Contact features count as 1.
 725    *
 726    * @param positional
 727    * @return
 728    */
 729
 730   public int getFeatureCount(boolean positional)
 731   {
 732     if (!positional)
 733     {
 734       return nonPositionalFeatures == null ? 0
 735               : nonPositionalFeatures.size();
 736     }
 737
 738     int size = 0;
 739
 740     if (contactFeatureStarts != null)
 741     {
 742       // note a contact feature (start/end) counts as one
 743       size += contactFeatureStarts.size();
 744     }
 745
 746     if (features != null)
 747     {
 748       size += features.size();
 749     }
 750
 751     return size;
 752   }
 753
 754   /**
 755    * Answers the total length of positional features (or zero if there are
 756    * none). Contact features contribute a value of 1 to the total.
 757    *
 758    * @return
 759    */
 760
 761   public int getTotalFeatureLength()
 762   {
 763     return totalExtent;
 764   }
 765
 766   /**
 767    * Answers the minimum score held for positional or non-positional features.
 768    * This may be Float.NaN if there are no features, are none has a non-NaN
 769    * score.
 770    *
 771    * @param positional
 772    * @return
 773    */
 774
 775   public float getMinimumScore(boolean positional)
 776   {
 777     return positional ? positionalMinScore : nonPositionalMinScore;
 778   }
 779
 780   /**
 781    * Answers the maximum score held for positional or non-positional features.
 782    * This may be Float.NaN if there are no features, are none has a non-NaN
 783    * score.
 784    *
 785    * @param positional
 786    * @return
 787    */
 788
 789   public float getMaximumScore(boolean positional)
 790   {
 791     return positional ? positionalMaxScore : nonPositionalMaxScore;
 792   }
 793
 794   /**
 795    * Answers a list of all either positional or non-positional features whose
 796    * feature group matches the given group (which may be null)
 797    *
 798    * @param positional
 799    * @param group
 800    * @return
 801    */
 802
 803   public List<SequenceFeature> getFeaturesForGroup(boolean positional,
 804           String group)
 805   {
 806     List<SequenceFeature> result = new ArrayList<>();
 807
 808     /*
 809      * if we know features don't include the target group, no need
 810      * to inspect them for matches
 811      */
 812     if (positional && !positionalFeatureGroups.contains(group)
 813             || !positional && !nonPositionalFeatureGroups.contains(group))
 814     {
 815       return result;
 816     }
 817
 818     List<SequenceFeature> sfs = positional ? getPositionalFeatures()
 819             : getNonPositionalFeatures();
 820     for (SequenceFeature sf : sfs)
 821     {
 822       String featureGroup = sf.getFeatureGroup();
 823       if (group == null && featureGroup == null
 824               || group != null && group.equals(featureGroup))
 825       {
 826         result.add(sf);
 827       }
 828     }
 829     return result;
 830   }
 831
 832   /**
 833    * Adds the shift amount to the start and end of all positional features whose
 834    * start position is at or after fromPosition. Returns true if at least one
 835    * feature was shifted, else false.
 836    *
 837    * @param fromPosition
 838    * @param shiftBy
 839    * @return
 840    */
 841
 842   public synchronized boolean shiftFeatures(int fromPosition, int shiftBy)
 843   {
 844     /*
 845      * Because begin and end are final fields (to ensure the data store's
 846      * integrity), we have to delete each feature and re-add it as amended.
 847      * (Although a simple shift of all values would preserve data integrity!)
 848      */
 849     boolean modified = false;
 850     for (SequenceFeature sf : getPositionalFeatures())
 851     {
 852       if (sf.getBegin() >= fromPosition)
 853       {
 854         modified = true;
 855         int newBegin = sf.getBegin() + shiftBy;
 856         int newEnd = sf.getEnd() + shiftBy;
 857
 858         /*
 859          * sanity check: don't shift left of the first residue
 860          */
 861         if (newEnd > 0)
 862         {
 863           newBegin = Math.max(1, newBegin);
 864           SequenceFeature sf2 = new SequenceFeature(sf, newBegin, newEnd,
 865                   sf.getFeatureGroup(), sf.getScore());
 866           addFeature(sf2);
 867         }
 868         delete(sf);
 869       }
 870     }
 871     return modified;
 872   }
 873
 874   /**
 875    * Find all features containing this position.
 876    *
 877    * @param pos
 878    * @return list of SequenceFeatures
 879    * @author Bob Hanson 2019.07.30
 880    */
 881
 882   public List<SequenceFeature> findOverlappingFeatures(int pos,
 883           List<SequenceFeature> result)
 884   {
 885     if (result == null)
 886     {
 887       result = new ArrayList<>();
 888     }
 889
 890     if (contactFeatureStarts != null)
 891     {
 892       findContacts(contactFeatureStarts, pos, result, true);
 893       findContacts(contactFeatureEnds, pos, result, false);
 894     }
 895     if (featuresList != null)
 896     {
 897       findOverlaps(featuresList, pos, result);
 898     }
 899     return result;
 900   }
 901
 902   /**
 903    * Binary search for contact start or end at a given (Overview) position.
 904    *
 905    * @param l
 906    * @param pos
 907    * @param result
 908    * @param isStart
 909    *
 910    * @author Bob Hanson 2019.07.30
 911    */
 912   private static void findContacts(List<SequenceFeature> l, int pos,
 913           List<SequenceFeature> result, boolean isStart)
 914   {
 915     int low = 0;
 916     int high = l.size() - 1;
 917     while (low <= high)
 918     {
 919       int mid = (low + high) >>> 1;
 920       SequenceFeature f = l.get(mid);
 921       switch (Long.signum((isStart ? f.begin : f.end) - pos))
 922       {
 923       case -1:
 924         low = mid + 1;
 925         continue;
 926       case 1:
 927         high = mid - 1;
 928         continue;
 929       case 0:
 930         int m = mid;
 931         result.add(f);
 932         // could be "5" in 12345556788 ?
 933         while (++mid <= high && (f = l.get(mid)) != null
 934                 && (isStart ? f.begin : f.end) == pos)
 935         {
 936           result.add(f);
 937         }
 938         while (--m >= low && (f = l.get(m)) != null
 939                 && (isStart ? f.begin : f.end) == pos)
 940         {
 941           result.add(f);
 942         }
 943         return;
 944       }
 945     }
 946   }
 947
 948   BitSet bs = new BitSet();
 949
 950   /**
 951    * Double binary sort with bitset correlation
 952    *
 953    *
 954    * @param features
 955    * @param pos
 956    * @param result
 957    */
 958   private void findOverlaps(List<SequenceFeature> features, int pos,
 959           List<SequenceFeature> result)
 960   {
 961     int n = featuresList.size();
 962     if (n == 1)
 963     {
 964       checkOne(featuresList.get(0), pos, result);
 965       return;
 966     }
 967     if (orderedFeatureStarts == null)
 968     {
 969       rebuildArrays(n);
 970     }
 971     SequenceFeature sf = findClosestFeature(orderedFeatureStarts, pos);
 972     while (sf != null) {
 973       if (sf.end >= pos)
 974       {
 975         result.add(sf);
 976       }
 977       sf = sf.containedBy;
 978     }
 979   }
 980
 981   private void linkFeatures(SequenceFeature[] intervals)
 982   {
 983     if (intervals.length < 2)
 984     {
 985       return;
 986     }
 987     int maxEnd = intervals[0].end;
 988     for (int i = 1, n = intervals.length; i < n; i++)
 989     {
 990       SequenceFeature ithis = intervals[i];
 991       if (ithis.begin <= maxEnd)
 992       {
 993         ithis.containedBy = getContainedBy(intervals[i - 1], ithis);
 994       }
 995       if (ithis.end > maxEnd)
 996       {
 997         maxEnd = ithis.end;
 998       }
 999     }
1000   }
1001
1002   private SequenceFeature getContainedBy(SequenceFeature sf,
1003           SequenceFeature sf0)
1004   {
1005     int begin = sf0.begin;
1006     while (sf != null)
1007     {
1008       if (begin <= sf.end)
1009       {
1010         System.out.println("\nFS found " + sf0.index + ":" + sf0
1011                 + "\nFS in    " + sf.index + ":" + sf);
1012         return sf;
1013       }
1014       sf = sf.containedBy;
1015     }
1016     return null;
1017   }
1018
1019   private SequenceFeature findClosestFeature(SequenceFeature[] l, int pos)
1020   {
1021     int low = 0;
1022     int high = l.length - 1;
1023     while (low <= high)
1024     {
1025       int mid = (low + high) >>> 1;
1026       SequenceFeature f = l[mid];
1027       switch (Long.signum(f.begin - pos))
1028       {
1029       case -1:
1030         low = mid + 1;
1031         continue;
1032       case 1:
1033         high = mid - 1;
1034         continue;
1035       case 0:
1036
1037         while (++mid <= high && l[mid].begin == pos)
1038           {
1039             ;
1040           }
1041           mid--;
1042         return l[mid];
1043       }
1044     }
1045     // -1 here?
1046     return (high < 0 || low >= l.length ? null : l[high]);
1047   }
1048
1049   private void checkOne(SequenceFeature sf, int pos,
1050           List<SequenceFeature> result)
1051   {
1052     if (sf.begin <= pos && sf.end >= pos)
1053     {
1054       result.add(sf);
1055     }
1056     return;
1057   }
1058
1059   /*
1060    * contact features ordered by first contact position
1061    */
1062   private SequenceFeature[] orderedFeatureStarts;
1063
1064   private void rebuildArrays(int n)
1065   {
1066     if (startComp == null)
1067     {
1068       startComp = new StartComparator();
1069     }
1070     orderedFeatureStarts = new SequenceFeature[n];
1071
1072     for (int i = n; --i >= 0;)
1073     {
1074       SequenceFeature sf = featuresList.get(i);
1075       sf.index = i;
1076       orderedFeatureStarts[i] = sf;
1077     }
1078     Arrays.sort(orderedFeatureStarts, startComp);
1079     linkFeatures(orderedFeatureStarts);
1080   }
1081
1082   class StartComparator implements Comparator<SequenceFeature>
1083   {
1084
1085     int pos;
1086
1087     @Override
1088     public int compare(SequenceFeature o1, SequenceFeature o2)
1089     {
1090       int p1 = o1.begin;
1091       int p2 = o2.begin;
1092       return (p1 < p2 ? -1 : p1 > p2 ? 1 : 0);
1093     }
1094
1095   }
1096
1097   static StartComparator startComp;
1098
1099   // class EndComparator implements Comparator<SequenceFeature>
1100   // {
1101   //
1102   // int pos;
1103   //
1104   // @Override
1105   // public int compare(SequenceFeature o1, SequenceFeature o2)
1106   // {
1107   // int p1 = o1.end;
1108   // int p2 = o2.end;
1109   // int val = (p1 < p2 ? 1 : p1 > p2 ? -1 : 0);
1110   // return val;
1111   // }
1112   //
1113   // }
1114   //
1115   // static EndComparator endComp;
1116
1117 }