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