=head1 LICENSE

  Copyright (c) 1999-2012 The European Bioinformatics Institute and
  Genome Research Limited.  All rights reserved.

  This software is distributed under a modified Apache license.
  For license details, please see

    http://www.ensembl.org/info/about/code_licence.html

=head1 CONTACT

  Please email comments or questions to the public Ensembl
  developers list at <dev@ensembl.org>.

  Questions may also be sent to the Ensembl help desk at
  <helpdesk@ensembl.org>.

=cut

=head1 NAME

Bio::EnsEMBL::Utils::Collector

=head1 SYNOPSIS

  # Inherit this base module in your feature specific Collector
  # instance:

  package Bio::EnsEMBL::Funcgen::Collector::ResultFeature;
  use base('Bio::EnsEMBL::Utils::Collector');

  # ... and define package config variables
  $Bio::EnsEMBL::Funcgen::Collector::bin_model = 'SIMPLE';
  $Bio::EnsEMBL::Funcgen::Collector::window_sizes =
    [ 30, 65, 130, 260, 450, 648, 950, 1296 ];
  # Could replace 30 with 0 here for low density data at natural resolution

  $Bio::EnsEMBL::Utils::Collector::bin_method =
    'count';    # only used by collector

  $Bio::EnsEMBL::Utils::Collector::packed_size = 2;

  # ... or simply use this module in a script either defining package
  # config variables, or passing as parameters to the constructor:

  my $collector =
    Bio::EnsEMBL::Utils::BaseCollector->new( -pack_template => 'v' );

  $Bio::EnsEMBL::Funcgen::Collector::pack_template = 'v';

  # Config variables can also be over-ridden by passing a config hash to
  # the store_window_bins_by_Slice() method:

  $collector->store_window_bins_by_Slice( $slice, (
                                            -pack_template => 'v',
                                            -packed_size   => 2 ) );

  # NOTE: Over-riding default config variables can cause problems when
  # storing or fetching data. e.g. Fetch may revert to using defaults or
  # table partitions may not match window sizes.

=head1 DESCRIPTION

This package is the base Collector class which contains generic
getter/setter methods along with the main 'collecting' methods which
perform the majority of the work in generating compressed data
collections optimised for web display.  The bins produced are aimed at
storage in a BLOB representing an entire seq_region i.e. even bins with
no features/null data are encoded as a 0 score.  Non-BLOB collections
are currently not supported.

If your Collection class defines a Bio::EnsEMBL::Feature, then its
adaptor should inherit from the relevant Collection class.

The minimum prerequisites of the input features/data are that they have
a start() and end() method.  For instance a Bio::EnsEMBL::Features
generated from a database or parsed from a flat file.

NOTE: This Collector does not have a lightweight mode previously used
for dynamic/on the fly collecting i.e. it does not take advantage of
bypassing object creation via the related BaseFeatureAdaptor method.

=cut

package Bio::EnsEMBL::Utils::Collector;

use Bio::EnsEMBL::Utils::Argument  ('rearrange');
use Bio::EnsEMBL::Utils::Exception ('throw');
use strict;
use warnings;

### Global package config vars

# Defaults
our $max_view_width = 1000000;   # Max bp width in location/detailed view
our $max_data_type_size = 16777216;    # Default is 16MB for long BLOB
# This is really a guide value as this should be set in the inheriting
# Collector class by deducting the rest of the row size from this value.
# Is is upto the inheritor to handle checking whether this size has been 
# exceeded.

# NOTE: Theoretically the min window size is: slice_length/(16777216/2)
# So for human chr1:  249,250,621/(16777216/2) = 29.7 => 30.  However,
# this size does not seem to directly translate to the MySQL
# max_allowed_packet_size.  Increasing max_allowed_packet_size to 64MB
# solves this issue, and substr operation doesn't appear to incur any of
# the potential memory(4*) usage issues.

# Others global package variables which are set in the inheriting
# Collector class.
our ( $bin_model,   $bin_method, $pack_template,
      $packed_size, $window_sizes );

=head2 new

  Args       : None
  Example    :

    my $collector = Bio::EnsEMBL::XXX::Collector::FEATURE->new();
    $collector->store_windows_by_Slice($slice);

    # Where XXX is, e.g. Compara, FuncGen etc.

  Description: Simple new method to enable use of collector
               when not inherited by a descendant of
               Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor

  Returntype : Bio::EnsEMBL::XXX::Collector
  Exceptions : None
  Caller     : Collector script
  Status     : At Risk

=cut

sub new {
  return bless {}, $_[0];   # Simply blesses this class as an empty hash.

  # Do not set anything here, as will not be first in ISA for
  # FeatureAdaptors.  Hence, not guaranteed to be called.
}

=head2 new_assembly

  Args       : optional - string assembly version e.g. GRCh37
  Example    : $collector->new_assembly('GRCh37');
  Description: Getter/Setter for new assembly version which should be
               used to project only 0 wsize Collections.
  Returntype : string
  Exceptions : None
  Caller     : store_window_bins_by_Slice() or
               write_collection() in inheriting Collector class.
  Status     : At Risk

=cut

sub new_assembly {
  my ( $self, $new_asm ) = @_;

  if ( defined($new_asm) ) {
    $self->{'new_assembly'} = $new_asm;
  }

  return $self->{'new_assembly'};
}

### Setter/Getter methods for basic/mandatory config
#   Can also be set using package variables in the inheriting
#   Collector/Adaptor or run script.  Allows over-riding of defaults set
#   in Adaptor/Collector.

# Package variables used here instead of attrs to enable easy
# default config in inheriting class/script method. Provided
# for easy/standardised fetch access outside of this package
# i.e. Collectors/Adaptors

=head2 max_data_type_size

  Args       : optional - int  Maximum size of collection in bytes
  Example    : $collector->max_data_type_size($new_max_size);
  Description: Getter/Setter for max_data_type_size, default is
               currently set at in this class as 16777216 (16MB), for
               long BLOB.  This is used by the write_collection()
               method to determine when to build and store a compressed
               collection.
  Returntype : int
  Exceptions : None
  Caller     : bins_per_record() and
               write_collection() in inheriting Collector class.
  Status     : At Risk

=cut

sub max_data_type_size {
  my ( $self, $size ) = @_;

  # Validate is sensible integer

  if ( defined($size) ) {
    if ( $size !~ /^\d+$/ ) {
      throw("max_data_type_size must be a integer of bytes, not $size");
    }
    $max_data_type_size = $size;
  } elsif ( !defined($max_data_type_size) ) {
    # This should never happen as we have defaults in this module.
    throw(   'You must define a '
           . '$Bio::EnsEMBL::Utils::Collector::max_data_type_size '
           . 'or pass -max_data_type_size config' );
  }

  return $max_data_type_size;
}

=head2 max_view_width

  Args       : optional - int  Maximum width of view
  Example    : $collector->max_view_width($new_max_width);
  Description: Getter/Setter for max_view_width, default is currently
               set at in this class as 500000bp, for maximum level of
               zoom permitted by location view.
  Returntype : int
  Exceptions : None
  Caller     : general
  Status     : At Risk

=cut

sub max_view_width {
  my ( $self, $size ) = @_;

  # Validate is sensible integer

  if ( defined($size) ) {
    if ( $size !~ /^\d+$/ ) {
      throw("max_view_width must be a integer, not $size");
    }
    $max_view_width = $size;
  } elsif ( !defined $max_view_width ) {
    # This should never happen as we have defaults in this module.
    throw(   'You must define a '
           . '$Bio::EnsEMBL::Utils::Collector::max_view_width '
           . 'or pass -max_view_width config' );
  }

  return $max_view_width;
}

=head2 bin_method

  Args[0]    : optional - string name of bin method e.g. 'max_magnitude'
  Args[1]    : optional - Bio::EnsEMBL::Funcgen::Parsers::InputSet
  Example    : my $bin_method = $self->bin_method();
  Description: Getter/Setter for bin_method, default is normally set in
               the inheriting Collector class either by package variable
               or by passing a config hash via the store methods.
  Returntype : string
  Exceptions : Throws if cannot set by package variable
  Caller     : general
  Status     : At Risk

=cut

sub bin_method {
  my ( $self, $bmethod, $config ) = @_;

  if ( defined($bmethod) ) {
    $bin_method = $bmethod;
  }

  if ( !defined($bin_method) ) {
    throw(   'You must define a '
           . '$Bio::EnsEMBL::Utils::Collector::bin_method '
           . 'or pass -bin_method config' );
  }

  if ( !$self->can( "_calculate_" . $bin_method ) ) {
    throw("$bin_method is not a valid/available binning method");
  }

  my $set_up_method = "_set_up_" . $bin_method;
  if ( $self->can($set_up_method) ) {
    $self->$set_up_method($config);
  }

  return $bin_method;
}

=head2 bin_model

  Args       : optional - string bin model e.g. SIMPLE or COMPLEX
  Example    : my $bin_model = $self->bin_model;
  Description: Getter/Setter for bin_model, default should be set in
               inheriting Collector class.  Currently only supports
               'SIMPLE' bin model.
  Returntype : string
  Exceptions : Throws if bin_model is not SIMPLE
  Caller     : general
  Status     : At Risk

=cut

sub bin_model {
  my ( $self, $bmodel ) = @_;

  if ( defined($bmodel) ) {
    $bin_model = $bmodel;
  }

  if ( !defined($bin_model) ) {
    throw(   'You must define a '
           . '$Bio::EnsEMBL::Utils::Collector::bin_model '
           . 'or pass -bin_model config' );
  }

  if ( $bin_model ne 'SIMPLE' ) {
    throw(   'Bio::EnsEMBL::Utils::Collector does not yet support '
           . 'non-SIMPLE bin models' );
  }

  return $bin_model;
}


=head2 window_sizes

  Args       : optional - arrayref of window sizes
  Example    :

    foreach my $wsize ( @{ $collector->window_sizes } )
    {    # Do some collecting
    }

  Description: Getter/Setter for window_sizes.  Default should be set
               in inheriting Collector (if the config is dynamic),
               FeatureAdaptor class or script using package variable or
               this method.
               NOTE: Redefining these may cause a mismatch with the
               table partition definition.
  Returntype : arrayref of ints
  Exceptions : Throws if cannot set a valid array of int window sizes
  Caller     : general
  Status     : At Risk - rename bin_sizes?

=cut

sub window_sizes {
  my ( $self, $sizes ) = @_;

  if ( defined($sizes) ) {
    $window_sizes = $sizes;
  }

  if ( !(    ref($window_sizes)
          && ( ref($window_sizes) eq 'ARRAY' )
          && ( scalar(@$window_sizes) > 0 ) ) )
  {
    throw('Must pass -windows_sizes in the config '
        . 'or define $Bio::EnsEMBL::Utils::Collector::window_sizes '
        . 'in your Collector as an array ref of integer window_sizes' );
  }

  return $window_sizes;
}




=head2 has_window_size

  Args       : int - window size to validate
  Example    : if( $collector->has_window_size('30') ){
                   #Do something wrt to 30bp window size
               }

  Description: Simple utility method to validate whether this Collector
               has a given window_size
  Returntype : Boolean
  Exceptions : Throws if window size not specified
  Caller     : general
  Status     : At Risk

=cut


sub has_window_size{
  my ( $self, $size ) = @_;

  if(! defined $size){
	throw('You must pass a window size to validate');
  }

  return grep(/$size/, @$window_sizes); 
}




### Getter/Setters for BLOB collection config
# NOTE: Overriding the defaults here may cause a mismatch when the data
# is retrieved.

=head2 pack_template

  Args       : optional - string perl 'pack' template 
  Example    : $self->pack_template('v');
  Description: Getter/Setter for pack_template. Default should be set
               in inheriting Collector (if the config is dynamic),
               FeatureAdaptor class or script using package variable or
               this method.
  Returntype : string
  Exceptions : Throws if cannot set pack_template from package variable
  Caller     : FeatureAdaptor::_obj_from_sth
  Status     : At Risk

=cut

sub pack_template {
  my ( $self, $template ) = @_;

  if ( defined($template) ) {
    $pack_template = $template;
  }

  if ( !defined($pack_template) ) {
    throw(   'Must pass a per score '
           . '-pack_template in the config '
           . 'or define $Bio::EnsEMBL::Utils::Collector::pack_template '
           . 'in your Collector' );
  }

  return $pack_template;
}

=head2 packed_size

  Args       : optional - int size of perl 'pack' template in bytes
  Example    : $self->packed_size(2);
  Description: Getter/Setter for packed_size.  Default should be set
               in inheriting Collector (if the config is dynamic),
               FeatureAdaptor class or script using package variable or
               this method.
  Returntype : string
  Exceptions : Throws if cannot set pack_template from pacakge variable
  Caller     : current_packed_size() and
               FeatureAdaptor::_obj_from_sth()
  Status     : At Risk

=cut

sub packed_size {
  my ( $self, $size ) = @_;

  if ( defined($size) ) {
    $packed_size = $size;
  }

  if ( !defined($packed_size) ) {
    throw(   'Must pass -packed_size(wrt to pack_template) config '
           . 'or define $Bio::EnsEMBL::Utils::Collector::packed_size '
           . 'in your Collector' );
  }

  if ( $packed_size !~ /^\d+$/ ) {
    throw(   "$packed_size is not an integer, "
           . "must pass a size integer for packed_size "
           . "which specifies size of pack_template:\t"
           . $pack_template );
  }

  return $packed_size;
}

=head2 bins_per_record

  Example    : my $bin_per_records = $self->bin_per_record
  Description: Simple method to calculate the max number of bins
               allowed per record given the current config.
  Returntype : int
  Exceptions : None
  Caller     :
  Status     : At Risk

=cut

sub bins_per_record {
  return int( $max_data_type_size/$packed_size );
}


=head2 current_packed_size

  Arg[0]     : int - window size
  Example    : my $cps = $self->current_packed_size($wsize);
  Description: Simple method to calculate the max number of bins
               allowed per record given the current config.
  Returntype : int
  Exceptions : None
  Caller     :
  Status     : At Risk

=cut

sub current_packed_size {
  my ( $self, $wsize ) = @_;
  return ( scalar( @{ $self->score_cache($wsize) } )*$packed_size );
}


=head2 score_cache

  Arg[0]     : int - window size
  Example    : my $cps = $self->current_packed_size($wsize);
  Description: Handles caching of bin scores for each window size
  Returntype : arrayref
  Exceptions : Throws if no window size defined
  Caller     : current_packed_size() and store_collection()
  methods
  Status     : At Risk

=cut

sub score_cache {
  my ( $self, $wsize, $scores ) = @_;

  if ( !defined($wsize) ) {
    throw('Must pass a window size argument');
  }

  $self->{'score_cache'}{$wsize} ||= [];

  if ( defined($scores) ) {
    push( @{ $self->{'score_cache'}{$wsize} }, @{$scores} );
  }

  return $self->{'score_cache'}{$wsize};
}


=head2 collection_start

  Arg[0]     : int - window_size
  Arg[1]     : optional int - seq_region_start
  Example    : my $coll_start->(150);
  Description: Getter/Setter collection seq_region_start
  Returntype : int
  Exceptions : Throws if no window size defined
  Caller     : store_window_bin_by_Slice() and write_collection()
  Status     : At Risk

=cut

sub collection_start {
  my ( $self, $wsize, $sr_start ) = @_;

  if ( !defined($wsize) ) {
    throw('Must pass a window size argument');
  }

  if ( defined($sr_start) ) {
    $self->{'collection_start'}{$wsize} = $sr_start;
  }

  return $self->{'collection_start'}{$wsize};
}


=head2 collection_end

  Arg[0]     : int - window_size
  Arg[1]     : optional int - seq_region_end
  Example    : my $coll_end->(150);
  Description: Getter/Setter collection seq_region_end
  Returntype : int
  Exceptions : Throws if no window size defined
  Caller     : inheriting Collector write_collection method
  Status     : At Risk

=cut

sub collection_end{
  my ($self, $wsize, $sr_end) = @_;
  throw('Must pass a window size argument') if ! defined $wsize;

  if(defined $sr_end){
	$self->{'collection_end'}{$wsize} = $sr_end;
  }
  else{
	return $self->{'collection_end'}{$wsize};
  }
}


=head2 collection_strand

  Arg[0]     : int - window_size
  Arg[1]     : optional int - seq_region_strand
  Example    : my $coll_start->(0);
  Description: Getter/Setter collection seq_region_strand
  Returntype : int
  Exceptions : Throws if no window size defined
  Caller     : inheriting Collector write_collection method
  Status     : At Risk - Collections are currently strandless

=cut

sub collection_strand {
  my ( $self, $wsize, $strand ) = @_;

  if ( !defined($wsize) ) {
    throw('Must pass a window size argument');
  }

  if ( defined $strand ) {
    $self->{'collection_strand'}{$wsize} = $strand;
  }

  return $self->{'collection_strand'}{$wsize};
}


### Here follows the actual working methods

=head2 _get_Slice_chunks

  Description: Defines the optimal set of slice chunks to use for
               generating collections such that redundant fetches
               are minimized.
  Returntype : hashref of window_size chunk size pairs
  Exceptions : Throws if no window sizes or max_view_width defined
  Caller     : store_window_bin_by_Slice()
  Status     : At Risk

=cut

sub _get_Slice_chunks {
  my $self = shift;

  if ( !defined($window_sizes) || !defined($max_view_width) ) {
    throw(   'You must pass both a window_size array ref '
           . 'and max_view_width arguments' );
  }

  if ( !defined( $self->{'_slice_chunks'} ) ) {
    # Calulate sensible slice length based on window sizes
    my @wsizes = sort { $a <=> $b } @$window_sizes;

    # Handle calculating only 0 wsize
    if ( scalar(@wsizes) == 1
         && $wsizes[0] == 0 )
    {
      return { $max_view_width => [0] };
    }

    my $multiplier = int( $max_view_width/$wsizes[$#wsizes] );

    my $chunk_length  = $multiplier*$wsizes[$#wsizes];
    my $not_divisible = 1;

    my %chunk_windows;   # Registry of chunk lengths to run with windows
    my %workable_chunks = map { $_ => {} } @wsizes;

    # get rid of natural resolution as this will always work
    delete $workable_chunks{'0'};

    while ( $not_divisible && $chunk_length != 0 ) {
      $not_divisible = 0;

      foreach my $wsize (@wsizes) {
        if ( $wsize == 0 ) {
          # Special wsize for normal data
          next;
        }

        # Set not divisible if modulus is true
        if ( $chunk_length % $wsize ) {
          $not_divisible = 1;
        } else {
          $workable_chunks{$wsize}{$chunk_length} = [];
        }
      }

      # Gradually shrink the length until we find a workable slice
      # length for all windows.
      if ($not_divisible) {
        $chunk_length -= $wsizes[$#wsizes];
      }
    }

    my %chunk_sets;

    if ( $chunk_length == 0 ) {
      print "Could not find chunk length "
        . "for all window sizes, "
        . "attempting to subset windows "
        . "using alternate slice length\n";

      foreach my $wsize ( keys(%workable_chunks) ) {
        # Loop through windows, seeing if they are workable in the other
        # windows.

        foreach my $chunk ( keys( %{ $workable_chunks{$wsize} } ) ) {

          foreach my $other_wsize ( keys %workable_chunks ) {
            next if $wsize == $other_wsize;

            if ( exists( $workable_chunks{$other_wsize}{$chunk} ) ) {
              # only push it onto the other wsize, as we will do the
              # reverse later
              $chunk_sets{$chunk}{$wsize} = undef;
            }
          }
        }
      }

      # %chunk_sets represents co-occurence of wsizes with repect to
      # chunks.  Take the set which has the most windows and the longest
      # chunk.  Then get the largest which handles the rest.

      # define possible set lengths
      my $i = 0;
      my %set_lengths;

      map { $set_lengths{$i} = []; $i++ } @wsizes;

      # get rid of natural resolution as this will always work
      delete $set_lengths{'0'};

      # Store chunks lengths for each set size
      foreach my $chunk ( keys(%chunk_sets) ) {
        my $set_size = scalar( values( %{ $chunk_sets{$chunk} } ) );
        push( @{ $set_lengths{$set_size} }, $chunk );
      }

      # Get the biggest set with the longest length;

      # Scalar here as we are disregarding natural resolution of 0 in
      # loop.
      my $largest_size      = scalar(@wsizes);
      my $found_largest_set = 0;

      while ( !$found_largest_set ) {
        $largest_size--;

        if ( scalar( @{ $set_lengths{$largest_size} } ) > 0 ) {
          $found_largest_set = 1;
        }
      }

      my ($largest_chunk) =
        sort { $b <=> $a } @{ $set_lengths{$largest_size} };

      my @largest_windows = keys %{ $chunk_sets{$largest_chunk} };
      @{ $chunk_windows{$largest_chunk} } = @largest_windows;

      print "Largest chunk $largest_chunk($largest_size) "
        . "contains windows: @largest_windows\n";

      my %remaining_windows = map { $_ => {} } @wsizes;

      # get rid of natural resolution as this will always work
      delete $remaining_windows{'0'};

      map { delete $remaining_windows{$_} } @largest_windows;
      my $remaining_set_size = scalar( keys(%remaining_windows) );

      # Use array here for practicality, would need to maintain hash if
      # we need to iterate.
      my @rwindows = keys(%remaining_windows);

      # Could be one window, but this will not be in the co-occurence
      # hash %chunk_sets.
      my $next_chunk;

      if ( scalar(@rwindows) == 1 ) {
        my ($last_window) = @rwindows;
        # Find a suitably large chunk for this one window.
        $multiplier = int( 500000/$last_window );
        $next_chunk = $multiplier*$last_window;
      } else {

        foreach my $chunk ( sort { $b <=> $a }
                            @{ $set_lengths{$remaining_set_size} } )
        {
          my $seen_count = 0;

          foreach my $rwindow (@rwindows) {
            if ( grep /$rwindow/,
                 ( values( %{ $chunk_sets{$chunk} } ) ) )
            {
              $seen_count++;
            }
          }

          if ( $seen_count == $remaining_set_size ) {
            $next_chunk = $chunk;
            last;
          }
        }

      }

      @{ $chunk_windows{$next_chunk} } = @rwindows;

      if ( defined($next_chunk) ) {
        print "Found next chunk length $next_chunk "
          . "contains remaining windows:\t@rwindows\n";
      } else {
        warn "Need to write iterative method for set definition";
        throw(   'Could not find workable slice length '
               . 'for remaining windows: '
               . join( ', ', @rwindows ) );
      }
    } else {
      @{ $chunk_windows{$chunk_length} } = keys(%workable_chunks);
      print "Found workable chunk length $chunk_length  "
        . "for all window sizes:\t"
        . join( ' ', @{ $chunk_windows{$chunk_length} } ) . "\n";
    }

    $self->{'_slice_chunks'} = \%chunk_windows;
  } ## end if ( !defined( $self->...))

  return $self->{'_slice_chunks'};
} ## end sub _get_Slice_chunks




=head2 set_config

  Arg[0]     : optional hash - parameter hash(see above methods for more info):

        WINDOW_SIZES        => array ref - subset of defined window
                               sizes
        BIN_METHOD          => string
        MAX_VIEW_WIDTH      => int
        MAX_DATA_TYPE_SIZE  => int
        PACK_TEMPLATE       => string
        PACKED_SIZE         => int
        BIN_MODEL           => string
        NEW_ASSEMBLY        => string
        METHOD_CONFIG       => hash of method specific config params
        SKIP_ZERO_WINDOW    => boolean - skips generation of 0 wsize
                               this is used if already generated
                               from an assembly projection.

               NOTE: Over-riding any of the default config may cause
               problems when storing or retrieving Collection data,
               except sub sets of default window sizes.

  Description: This method replaces the constructor as new will not be
               called for Adaptor based Collectors.
               Separating this from the store method is currently 
               redundant as jobs are normally submitetd in Slice based 
               jobs. However, this will be required if the store method 
               is further seaprated into fetch/generate and store methods
  Returntype : None
  Exceptions : Throws if no window sizes or max_view_width defined
  Caller     : Inheritor Collector e.g. Bio::EnsEMBL::Funcgen:Collector::ResultFeature
               or script.
  Status     : At Risk

=cut

sub set_config {
  my ( $self, %config ) = @_;

  my ( $wsizes,       $bmethod,  $mv_width,
       $md_type_size, $template, $psize,
       $bmodel,       $new_assm, $skip_zero_window,
       $method_config )
    = rearrange( [ 'WINDOW_SIZES',     'BIN_METHOD',
                   'MAX_VIEW_WIDTH',   'MAX_DATA_TYPE_SIZE',
                   'PACK_TEMPLATE',    'PACKED_SIZE',
                   'BIN_MODEL',        'NEW_ASSEMBLY',
                   'SKIP_ZERO_WINDOW', 'METHOD_CONFIG' ],
                 %config );

  ### VAILDATE/SET VARS/CONFIG

  # Attrs used in this method
  $self->bin_method( $bmethod, $method_config );
  $self->bin_model($bmodel);
  $self->window_sizes($wsizes);

  # Set to undef if we have empty array? To change this we need to
  # pass the config hash -window_sizes conditionally
  # This currently overwrite the defaults!
  #  if ( ref($window_sizes) eq 'ARRAY'
  #       && scalar( @{$window_sizes} ) == 0 )
  #  {
  #    $window_sizes = undef;
  #  }

  # Attrs used in other (store) methods
  $self->pack_template($template);
  $self->packed_size($psize);
  $self->max_data_type_size($md_type_size);
  $self->max_view_width($mv_width);

  # Other vars
  $self->new_assembly($new_assm);
  $self->{'_only_natural'} = 0;
  $self->{'_store_natural'} = grep /^0$/, @$window_sizes;

  ### Set window_sizes

  if ( $self->new_assembly() ) {
    print "Assembly projection may cause problems "
      . "for large Collections, "
      . "defaulting to window_sizes = (0)\n";


	if ( $skip_zero_window ) {
	  throw(   "You cannot -skip_zero_window or "
			   . "omit 0 from -window_sizes "
			   . "when projecting to a new assembly($new_assm) "
			   . "which should only be generated using window_size=0" );
	}
	



    # Then build the bins on the projected 0 level single Features

    # Test we haven't explicity set window_sizes to be something else
    if ( defined($wsizes)
         && !( scalar(@$wsizes) == 1 && $wsizes->[0] == 0 ) )
    {
      throw(   "You have set window_sizes config "
             . "which are not safe when projecting to "
             . "a new assembly($new_assm), "
             . "please omit window_sizes config or set to 0" );
    }

    $self->window_sizes( [0] );
  } else {

    if ( $wsizes && $skip_zero_window && 
		 ( grep /^0$/, @$wsizes )) {
	  #Only test passed params not default config

      throw(   "You have specied skip_zero_window "
             . "and window_size 0 in your parameters, "
             . "please remove one of these" );
    } 
	elsif ( defined($window_sizes) && !grep /^0$/, @$window_sizes ) {
      $skip_zero_window = 1;
      # re-add 0 window as we need this to build the collections
	  # see ...
      unshift( @{$window_sizes}, 0 );
    }
  }

 
  if ( $self->{'_store_natural'} && scalar( @{$window_sizes} ) == 1 ) {
    $self->{'_only_natural'} = 1;
  }
  if ($skip_zero_window) {
    $self->{'_store_natural'} = 0;
  }

  return;
} ## end sub set_config

=head2 store_window_bins_by_Slice

  Arg[0]     : Bio::EnsEMBL:Slice
  Example    : $collector->store_window_bins_by_Slice($slice);
  Description: This is the main run method, it loops through
               optimal slice chunks from _define_window_chunks,
               calls _bin_features_by_Slice as appropriate and
               calls write_collection in the inheriting Collector
               class/script.
  Returntype : None
  Exceptions : Throws if Bio::EnsEMBL::Slice is not defined
  Caller     : store methods in inheriting Collector class/script
  Status     : At Risk

=cut

sub store_window_bins_by_Slice {
  my ( $self, $slice ) = @_;

  warn "Need to be careful here "
    . "about cleaning start end strand caches between "
    . "serially run slices";

  if ( !(    defined($slice)
          && ref($slice)
          && $slice->isa('Bio::EnsEMBL::Slice') ) )
  {
    throw('You must pass a valid Bio::EnsEMBL::Slice');
  }

  # Rollback previously stored features.
  # Change 'can' to empty method stubb with pod ???
  if ( $self->can('rollback_Features_by_Slice') ) {
    $self->rollback_Features_by_Slice($slice);
  } else {
    warn ref($self)
      . " cannot rollback_Features_by_Slice. "
      . "This may result in storage failure "
      . "or duplicate Collections if there is pre-existing data";
  }

  ### PROCESS CHUNKS
  my %chunk_windows = %{ $self->_get_Slice_chunks };
  my (%counts);
  my $store_natural = $self->{'_store_natural'};
  my $only_natural  = $self->{'_only_natural'};
  $counts{0} = 0;    # Set natural res count to 0
  my $slice_end       = $slice->end;
  my $orig_start      = $slice->start;
  my $region          = $slice->coord_system_name;
  my $version         = $slice->coord_system->version;
  my $seq_region_name = $slice->seq_region_name;
  my $strand          = $slice->strand;

  # Warn if this is not a full slice.  Version needed in case we are
  # projecting from a non-default version slice
  my $full_slice =
    $slice->adaptor->fetch_by_region( $region, $seq_region_name, undef,
                                      undef, undef, $version );

  if (    ( $full_slice->start() != $orig_start )
       || ( $full_slice->end() != $slice_end ) )
  {
    warn "Generating collections using sub-Slices "
      . "can result in data issues/artifacts";
    # Last chunk might not be the correct window length.  Test
    # slices less than chunk length can cause failures in
    # _bin_features_by_window_sizes others?
  }

  # Set the initial collection_start to orig_start.  This is not the
  # case for 0 wsize where it must always be the true feature start.
  for my $wsize (@$window_sizes) {
    if ( $wsize == 0 ) { next }
    $self->collection_start( $wsize, $orig_start );

    # Also reset collection end and score cache in case we are running
    # serially.
    $self->{collection_end}{$wsize} = undef;
    $self->{'score_cache'}{$wsize} = [];
  }

  my $first_chunk_length = 1;

  foreach my $chunk_length ( sort keys %chunk_windows ) {
    print "Processing windows "
      . join( ', ', @{ $chunk_windows{$chunk_length} } )
      . " with chunk length $chunk_length\n";

    # Set window counts to 0
    map $counts{$_} = 0, @{ $chunk_windows{$chunk_length} };

    # May need to reset flat file parser handle or other caches via
    # inheriting Collector
    if ( !$first_chunk_length ) {
      # Change 'can' to empty method stubb with pod???
      if ( $self->can('reinitialise_input') ) {
        $self->reinitialise_input();
      }
    }

    $first_chunk_length = 0;

    # Now walk through slice using slice length chunks and build all
    # windows in each chunk.
    my $in_slice  = 1;
    my $start_adj = 0;
    my ( $sub_slice, $sub_end, $features, $bins );
    my $sub_start    = 1;
    my $slice_length = $slice->length();

    # Always create in local coords for fetch
    # Then change to seq_region coords for store if required

    while ($in_slice) {
      $sub_start += $start_adj;
      $sub_end = $sub_start + $chunk_length - 1;

      if ( $sub_end >= $slice_length ) {
        # Surplus bins are removed in store/write_collection in caller
        $in_slice = 0;
      }

      $sub_slice =
        $slice->adaptor->fetch_by_region( $region, $seq_region_name,
                                          $sub_start + $orig_start - 1,
                                          $sub_end + $orig_start - 1,
                                          $strand, $version );

      # Can't subslice as this will not clip if we go over the length of
      # the slice, unlike normal slice fetching.  Will clipping the end
      # to the slice end cause any problems here?  How will this affect
      # bin clipping?

      ### Grab features and shift chunk coords
      $features = $self->get_Features_by_Slice($sub_slice);

      # warn "Binning "
      #   . scalar(@$features)
      #   . " Features for chunk length $chunk_length, on Slice "
      #   . $sub_slice->name;

      if ( ( @{$features} )
         && ref( $features->[0] ) =~ /Bio::EnsEMBL::Utils::Collection/ )
      {
        # Would need to create base module with generic methods:
        # window_size, ...

      # Check that the returned feature/collections support window_size.
      # All Collections should be able to

        if ( $features->[0]->can('window_size') ) {
          if ( $features->[0]->window_size != 0 ) {
            throw(   "You are trying to generated Collections from "
                   . "a non-zero window sized Collection:\t"
                   . $features->[1]->{'window_size'} );
          }

        # This should never happen
        # if ( !$skip_zero_window ) {
        #   throw( 'You have retrieved data from a Collection '
        #        . 'which without using -skip_zero_window '
        #        . 'i.e. you are trying to generate overwrite '
        #        . 'the data you are generating the Collections from' );
        # }

        } else {
          throw(   'Something is wrong, '
                 . 'the Collection you have retrieved '
                 . 'does not support the method window_size' );
        }
      } ## end if ( ( @{$features} ) ...)

      # Set collection start here for 0 window_size
      if (    @{$features}
           && $store_natural
           && !defined( $self->collection_start(0) ) )
      {
        $self->collection_start( 0,
                                 $features->[0]->start + $sub_start );
      }

      if ($in_slice) {
        $start_adj = $chunk_length;
      }

      # Collect features into wsize bins
      if ( !$only_natural ) {
        # Get hashref of wsize=>bin array pairs
        $bins =
          $self->_bin_features_by_Slice_window_sizes(
                         -slice        => $sub_slice,
                         -window_sizes => $chunk_windows{$chunk_length},
                         -features     => $features, );
      }

      # Handle 0 wsize
      if ($store_natural) {
        foreach my $feature ( @{$features} ) {
          $counts{0}++;

          if ( $bin_model eq 'SIMPLE' ) {
            $self->collection_start( 0, $feature->start + $sub_start );

            $self->write_collection(
              0,
              $slice,    # Pass Slice to sub-slice when storing
              $feature->end + $sub_start,
              $feature->strand,   # Need to pass strand for 0 resolution
              $feature->scores, );
          }
        }

        print "Window size 0 (natural resolution) has "
          . scalar( @{$features} )
          . " feature bins for:\t"
          . $sub_slice->name . "\n";
      }

      # Now store collections for wsizes >0
      my $num_bins;

      foreach my $wsize ( sort keys( %{$bins} ) ) {
        $num_bins = scalar( @{ $bins->{$wsize} } );
        $counts{$wsize} += $num_bins;

        if ( $bin_model eq 'SIMPLE' ) {
          $self->write_collection(
            $wsize,
            $slice,
            #$sub_start,
            $sub_end,
            $slice->strand,    # This is most likely 1!
             # Override this woth 0 in descendant Collector if required.
            $bins->{$wsize}, );

        } else {
          throw(   'Bio::EnsEMBL::Utils::Collector '
                 . 'does not yet support non-SIMPLE bin models' );
          # i.e. More than one score
        }
      }
    } ## end while ($in_slice)

    # Turn off storing of natural resolution for next chunk length sets
    $store_natural = 0;
  } ## end foreach my $chunk_length ( ...)

  # Write last collections for each wsize

  foreach my $wsize (@$window_sizes) {

    if (    ( $wsize == 0 && !$store_natural )
         || ( $wsize != 0 && $only_natural ) )
    {
      next;
    }

    print "Writing final $wsize window_size collection, "
      . "this may result in slightly different "
      . "bin numbers from counts due to removing "
      . "overhanging bins past end of slice\n";
    $self->write_collection( $wsize, $slice );
  }

  # Print some counts
  foreach my $wsize ( sort ( keys %counts ) ) {
    print "Generated "
      . $counts{$wsize}
      . " bins for window size $wsize for "
      . $slice->name . "\n";
    # Some may have failed to store if we are projecting to a new
    # assembly.
  }

  return;
} ## end sub store_window_bins_by_Slice

=head2 _bin_features_by_Slice_window_sizes

  Args[0]    : Bio::EnsEMBL::Slice
  Args[1]    : ARRAYREF of window sizes
  Args[2]    : ARRAYREF of features with start and end method
               e.g. Bio::EnsEMBL::Features
  Example    :

      $bins =
        $self->_bin_features_by_window_sizes(
                         -slice        => $slice,
                         -window_sizes => $chunk_windows{$chunk_length},
                         -features     => $features, );

  Description: Bins feature scores for a given list of window sizes and
               predefined method.
  Returntype : HASHREF of scores per bin per window size
  Exceptions : None
  Caller     : store_window_bins_by_Slice
  Status     : At Risk

=cut

sub _bin_features_by_Slice_window_sizes {
  my ( $self, @args ) = @_;

  my ( $slice, $wsizes, $features ) =
    rearrange( [ 'SLICE', 'WINDOW_SIZES', 'FEATURES' ], @args );

  # Generate these once in caller?
  my $calc_method = '_calculate_' . $bin_method;
  my $post_method = '_post_process_' . $bin_method;

  # Do this conditional on the Collection type i.e. is
  # collection seq_region blob then no else yes Would need
  # $Bio::EnsEMBL::Utils::Collector::collection_format=BLOB|STANDARD
  # if ( !defined($features) || !@{$features} ) { return {} }

  # Set up some hashes to store data by window_size
  my ( %bins, %nbins, %bin_counts );
  my $slice_start  = $slice->start();
  my $slice_length = $slice->length();

  # Set up some bin data for the windows
  foreach my $wsize (@$wsizes) {
    $nbins{$wsize} = int( $slice_length/$wsize );    # int rounds down
         # nbins is index of the bin not the 'number'
         # Unless $slice_length is a multiple!
    if ( !( $slice_length % $wsize ) ) { $nbins{$wsize}-- }

    # Create default bins with 0
    $bins{$wsize} = [];
    map { $bins{$wsize}->[$_] = 0 } ( 0 .. $nbins{$wsize} );

    # Set bin counts to 0 for each bin
    $bin_counts{$wsize} = [];

    # This is adding an undef to the start of the array!?
    map { $bin_counts{$wsize}->[ ($_) ] = 0 } @{ $bins{$wsize} };

    foreach my $bin ( @{ $bins{$wsize} } ) {
      $bin_counts{$wsize}->[$bin] = 0;
    }
  }

  my $feature_index = 0;
  my ( $bin_index, @bin_masks );

  foreach my $feature ( @{$features} ) {
    # Set up the bins for each window size

    foreach my $wsize (@$wsizes) {
      my $start_bin = int( ( $feature->start )/$wsize );
      my $end_bin   = int( ( $feature->end )/$wsize );

      if ( $end_bin > $nbins{$wsize} ) {
        $end_bin = $nbins{$wsize};
      }

      $self->$calc_method( $feature, $start_bin, $end_bin,
                           $wsize,   \%bins,     \%bin_counts );
    }

  }

  # Now do post processing of bins if required
  if ( $self->can($post_method) ) {
    $self->$post_method( \%bins, \%bin_counts );
  }

  return \%bins;
} ## end sub _bin_features_by_Slice_window_sizes
# end sub _bin_features_by_Slice


### Here follows the bin methods
# These may also be defined in the inheriting Collector class.  No tests
# as these are internal and require speed.


=head2 _calculate_count

  Args[0]    : feature e.g. Bio::EnsEMBL::Feature
  Args[1]    : int - start bin
  Args[2]    : int - end bin
  Args[3]    : int - window_size
  Args[4]    : hashref - score bins
  Example    : $self->$calc_method
  Description: Adds count to bins which this feature overlaps
  Returntype : None
  Exceptions : None
  Caller     : _bin_features_by_window_sizes
  Status     : At Risk

=cut

sub _calculate_count {
  my ( $self, $feature, $start_bin, $end_bin, $wsize, $bins_ref ) = @_;

  my $bin_index;

  for ( $bin_index = $start_bin; $bin_index <= $end_bin; ++$bin_index )
  {
    $bins_ref->{$wsize}->[$bin_index]++;
  }

  return;
}


=head2 _calculate_average_score

  Args[0]    : feature e.g. Bio::EnsEMBL::Feature
  Args[1]    : int - start bin
  Args[2]    : int - end bin
  Args[3]    : int - window_size
  Args[4]    : hashref - score bins
  Example    : $self->$calc_method
  Description: Adds score to bins which this feature overlaps
  Returntype : None
  Exceptions : None
  Caller     : _bin_features_by_window_sizes
  Status     : At Risk

=cut


sub _calculate_average_score {
  my ( $self, $feature, $start_bin, $end_bin, $wsize, $bins_ref,
       $bin_counts_ref )
    = @_;

  # This is simple an average of all the scores for features which
  # overlap this bin.  No weighting with respect to the bin or the
  # feature.

  my $score = $self->get_score_by_Feature($feature);

  for ( my $bin_index = $start_bin;
        $bin_index <= $end_bin;
        ++$bin_index )
  {
    # We should really push onto array here so we can have median or
    # mean.

    $bins_ref->{$wsize}->[$bin_index] += $score;
    $bin_counts_ref->{$wsize}->[$bin_index]++;
  }

  return;
}


=head2 _post_process_average_score

  Args[0]    : hashref - score bins
  Args[1]    : hashref - count bins
  Example    : $self->$post_method
  Description: Post processes bins to calculate average score
  Returntype : None
  Exceptions : None
  Caller     : _bin_features_by_window_sizes
  Status     : At Risk

=cut	

sub _post_process_average_score {
  my ( $self, $bins_ref, $bin_counts_ref ) = @_;

  foreach my $wsize ( keys %{$bins_ref} ) {
    foreach my $bin_index ( 0 .. $#{ $bins_ref->{$wsize} } ) {

      if ( $bin_counts_ref->{$wsize}->[$bin_index] ) {
        $bins_ref->{$wsize}->[$bin_index] /=
          $bin_counts_ref->{$wsize}->[$bin_index];
      }

    }
  }

  return;
}


=head2 _calculate_max_magnitude

  Args[0]    : feature e.g. Bio::EnsEMBL::Feature
  Args[1]    : int - start bin
  Args[2]    : int - end bin
  Args[3]    : int - window_size
  Args[4]    : hashref - score bins
  Example    : $self->$calc_method
  Description: Sets max +/-ve scores for bins which this feature overlaps
  Returntype : None
  Exceptions : None
  Caller     : _bin_features_by_window_sizes
  Status     : At Risk

=cut

sub _calculate_max_magnitude {
  my ( $self, $feature, $start_bin, $end_bin, $wsize, $bins_ref ) = @_;

  my $score = $self->get_score_by_Feature($feature);

  # Max magnitude
  # Take the highest value +ve or -ve score
  for ( my $bin_index = $start_bin;
        $bin_index <= $end_bin;
        ++$bin_index )
  {

    # We really need to capture the lowest -ve and higest +ve scores
    # here and post process to pick between them.

    $bins_ref->{$wsize}->[$bin_index] ||= [ 0, 0 ];    #-ve, +ve

    if ( $score < $bins_ref->{$wsize}->[$bin_index]->[0] ) {
      $bins_ref->{$wsize}->[$bin_index]->[0] = $score;
    } elsif ( $score > $bins_ref->{$wsize}->[$bin_index][1] ) {
      $bins_ref->{$wsize}->[$bin_index]->[1] = $score;
    }
  }

  return;
} ## end sub _calculate_max_magnitude


=head2 _post_process_max_magnitude

  Args[0]    : hashref - score bins
  Args[1]    : hashref - count bins
  Example    : $self->$post_method
  Description: Post processes bins to pick largest +ve or -ve score
  Returntype : None
  Exceptions : None
  Caller     : _bin_features_by_window_sizes
  Status     : At Risk

=cut	

sub _post_process_max_magnitude {
  my ( $self, $bins_ref ) = @_;

  # Take the highest value +ve or -ve score

  foreach my $wsize ( keys %{$bins_ref} ) {
    foreach my $bin_index ( 0 .. $#{ $bins_ref->{$wsize} } ) {

      # Have potential for no listref in a given bin

      # default value if we haven't seen anything is 0
      # Actually want an array of -ve +ve values

      if ( $bins_ref->{$wsize}->[$bin_index] ) {
        my $tmp_minus = -$bins_ref->{$wsize}->[$bin_index]->[0];

        if ( $tmp_minus > $bins_ref->{$wsize}->[$bin_index]->[1] ) {
          $bins_ref->{$wsize}->[$bin_index] =
            $bins_ref->{$wsize}->[$bin_index]->[0];
        } else {
          $bins_ref->{$wsize}->[$bin_index] =
            $bins_ref->{$wsize}->[$bin_index]->[1];
        }

      }

    }
  }
  return;
} ## end sub _post_process_max_magnitude


=head2 _calculate_RPKM

  Args[0]    : feature e.g. Bio::EnsEMBL::Feature
  Args[1]    : int - start bin
  Args[2]    : int - end bin
  Args[3]    : int - window_size
  Args[4]    : hashref - score bins
  Example    : $self->$calc_method
  Description: Stores counts to calculate Read Per Kb per Million(RPKM)
  Returntype : None
  Exceptions : None
  Caller     : _bin_features_by_window_sizes
  Status     : At Risk

=cut

sub _calculate_RPKM {
  my ( $self, $feature, $start_bin, $end_bin, $wsize, $bins_ref ) = @_;

  $self->_calculate_count( $feature, $start_bin, $end_bin,
                           $wsize,   $bins_ref );

  return;
}


=head2 _post_process_RPKM

  Args[0]    : hashref - score bins
  Args[1]    : hashref - count bins
  Example    : $self->$post_method
  Description: Post processes bins to calculate average score
  Returntype : None
  Exceptions : None
  Caller     : _bin_features_by_window_sizes
  Status     : At Risk

=cut	

sub _post_process_RPKM {
  my ( $self, $bins_ref ) = @_;

  #10^9 x C / NGB
  #C = Reads overlapping bin
  #N = Total reads in the experiment
  #G = Length of bin in bps 
  #(don't really have to account for non-ref/HAPs or gender here 
  #as should be close enough, CellTypes/gender differences will be miniscule)
  #B = length of each bin

  foreach my $wsize(keys %{$bins_ref}){
	
	foreach my $bin_index(0..$#{$bins_ref->{$wsize}}){
	  $bins_ref->{$wsize}->[$bin_index] =  
		((10**9) *  
		 $bins_ref->{$wsize}->[$bin_index])/(($self->_RPKM_factor($wsize)) * $wsize);
	}
  }

  return;

}


=head2 _set_up_RPKM

  Args[0]    : hashref - method config e.g
                 {
                  DNADB         => Bio::EnsEMBL::DBSQL::DBAdaptor,
                  TOTAL_FEATURE => $total_feature_count,
                 }

  Example    : $self->$set_up_method($config);
  Description: Sets the RPKM factor
  Returntype : None
  Exceptions : Throws is required config params are not set
  Caller     : bin_method
  Status     : At Risk

=cut	


sub _set_up_RPKM{
  my ($self, $config) = @_;

  my ($dnadb, $total_features) = rearrange([ 'DNADB', 'TOTAL_FEATURES'], %{$config});
  
  #Test specifically here to notify about config hash
  if(! $total_features){
	throw("For RPKM you must pass a valid 'total_features' ".
		  "as part of the method config hash.");
  }
  
  if(! $dnadb){
	throw("For RPKM you must pass 'dnadb' as part of the method config hash.");
  }

  foreach my $wsize(@{$self->window_sizes}){
	#Should never have 0 here
	$self->_RPKM_factor($wsize, ($wsize * $total_features));	#N*G

	warn "setting $wsize RPKM factor($wsize * $total_features) to ".	
	  $self->_RPKM_factor($wsize);
  }
  
  return;
} ## end sub _set_up_RPKM


=head2 _RPKM_factor

  Args[0]    : int - RPKM factor i.e.  (Total reads in the experiment *
               Genome length)
  Example    : $self->_RPKM_factor($wsize, $factor);
  Description: Gets/Sets the RPKM factor
  Returntype : int
  Exceptions : None
  Caller     : _set_up_RPKM, _post_process_RPKM
  Status     : At Risk

=cut	

sub _RPKM_factor{
  my ($self, $wsize, $factor) = @_;

  if (! defined $wsize){
	throw("You must pass at least window_size to get or set the RPKM factor");
  }

  if(defined $factor){
	$self->{'RPKM_factor'}{$wsize} = $factor;
  }
  elsif(! exists $self->{'RPKM_factor'}{$wsize}){
	#This should never happen unless the window sizes 
	#are redefined after initialisation
	throw("You have requested an RPKM factor for a window_size".
		  " which has not been set:\t$wsize");
  }

  return $self->{'RPKM_factor'}{$wsize};
}

=head2 get_diploid_genome_length_by_gender

  Args[0]    : string - RPKM factor i.e.  (Total reads in the experiment *
               Genome length)
  Args[1]    : string - gender e.g. male or female
  Example    :

      my $glength =
        $self->get_diploid_genome_length_by_gender( $dnadb, $gender );

  Description: Gets the gender specific diploid genome length,
               including non-ref but not including haplotypes.  Only
               handles species with X/Y sex chromosomes.
  Returntype : int
  Exceptions : None
  Caller     : _set_up_RPKM, _post_process_RPKM
  Status     : At Risk - Move to and export from generic Utils Slice module???

=cut	

sub get_diploid_genome_length_by_gender {
  my ( $dnadb, $gender ) = @_;

  my %sex_chrs = ( 'Y' => 'male',
                   'X' => 'female', );

  my $dip_length = 0;

  if (!(( ref($dnadb) && $dnadb->isa('Bio::EnsEMBL::DBSQL::DBAdaptor') )
        && $dnadb->grou() eq 'core'
        && ( defined $gender && $gender =~ /(male|female)/ ) ) )
  {
    throw(   "Must provide valid "
           . "Bio::EnsEMBL::DBSQL::DBAdaptor($dnadb) and "
           . "gender ($gender) arguments" );
  }

  my @ref_slices = $dnadb->get_SliceAdaptor->fetch_all('toplevel');

  # Include non-ref(unassembled), but omit haps/lrgs(i.e. redundant)

  foreach my $slice (
     @{ $dnadb->get_SliceAdaptor->fetch_all( 'toplevel', undef, 1, 1 ) }
    )
  {
    # Include duplicated region for true diploid length

    # Skip haps/lrgs
    if ( ( $slice->coord_system->name() eq 'chromosome'
           && !$slice->is_reference() )
         || $slice->coord_system->name() eq 'lrg' )
    {
      next;
    }

    if ( exists( $sex_chrs{ $slice->seq_region_name() } ) ) {
      if ( $gender eq 'male' ) {
        $dip_length += $slice->length;
      } elsif ( $sex_chrs{ $slice->seq_region_name } eq 'male' ) {
        next;
      }
    }

    $dip_length += 2*$slice->length;
  }

  return $dip_length;
} ## end sub get_diploid_genome_length_by_gender


1;
