# $Id: ChaosGraph.pm,v 1.7 2005/06/15 16:21:09 cmungall Exp $
#
#

=head1 NAME

  Bio::Chaos::ChaosGraph     - object for representing a chaos-xml dataset

=head1 SYNOPSIS

  use Bio::Chaos::ChaosGraph;
  use Data::Stag;

  my $chaos_node = Data::Stag->parse('Rab1.chaos');
  my $cg = Bio::Chaos::ChaosGraph->new($chaos_node);
  
  my $fl = $cg->top_features;
  foreach my $f (@$fl) {
    next unless $f->get_type eq 'gene';
    $island_feature = $cg->make_gene_island($f, 5000, 5000);

    print $island_feature->xml;
  }

=head1 DESCRIPTION

This class allows manipulation of in-memory Chaos documents as
L<Data::Stag> nodes, and provides additional methods for traversing
the graph structure defined in the Chaos document

=head1 SEE ALSO

The L<cx-genbank2chaos.pl> script

The L<Bio::Chaos> helper module

The BioPerl L<Bio::SeqIO::chaosxml> writer

=cut

package Bio::Chaos::ChaosGraph;

use Exporter;
use Data::Stag qw(:all);
use Bio::Chaos::Root;
@ISA = qw(Bio::Chaos::Root Exporter);

use FileHandle;
use strict;
use Graph;

# Constructor


=head2 new

  Usage   - my $chaos = Bio::Chaos::ChaosGraph->new($chaos_stag)
  Returns - Bio::Chaos::ChaosGraph

creates a new Chaos::ChaosGraph object

=cut

sub new {
    my $proto = shift; my $class = ref($proto) || $proto;;
    my $self = bless {}, $class;
    my ($stag,$file,$fmt) =
      $self->_rearrange([qw(stag file fmt)], @_);
    if ($stag && $file) {
        $self->freak("use -stag OR -file as arguments");
    }
    #my @g_opt = (compat02=>1);
    my @g_opt = ();
    $self->graph(Graph->new(@g_opt));
    $self->locgraph(Graph->new(@g_opt));
    $self->feature_idx({});
    $self->init_from_stag($stag) if $stag;
    $self->init_from_file($file,$fmt) if $file;
    return $self;
}


=head2 init_from_stag

  Usage   - $cg->init_from_stag($chaos_node);
  Returns -
  Args    - L<Data::Stag>

sets up a CG object from a stag node conforming to the Chaos-XML schema

  my $cg = Bio::Chaos::ChaosGraph->new;
  my $cn = Data::Stag->parse("mydata.chaos.xml");
  $cg->init_from_stag($cn);

=cut

sub init_from_stag {
    my $self = shift;
    my $stag = shift;
    if (!$stag) {
	$self->freak;
    }
    my $verbose = $self->verbose;
    $self->debug("Adding features") if $verbose; 
    foreach my $feature ($stag->get_feature) {
	$self->add_feature($feature);
    }
    $self->debug("Features added") if $verbose; 
    $self->debug("Adding feature_relationships") if $verbose; 
    foreach my $fr ($stag->get_feature_relationship) {
	$self->add_feature_relationship($fr);
    }
    $self->debug("Feature_relationships added") if $verbose; 
    return;
}

=head2 init_from_file

  Usage   - $cg->init_from_file($chaos_node);
  Returns -
  Args    - filename string

sets up a CG object from a file conforming to the Chaos-XML schema

  my $cg = Bio::Chaos::ChaosGraph->new;
  $cg->init_from_file("mydata.chaos.xml");

=cut

sub init_from_file {
    my $self = shift;
    my $file = shift;
    my $fmt = shift;
    if (!$fmt) {
        if ($file =~ /chaos/) {
            $fmt = 'chaos';
        }
        else {
            $fmt = 'genbank';
        }
    }
    if ($fmt eq 'chaos') {
        my $stag = Data::Stag->parse($file);
	$self->init_from_stag($stag);
        return;
    }
    $self->chaos_flavour("$fmt-unflattened");
    $self->load_module("Bio::SeqIO");
    my $unflattener = $self->unflattener;
    my $type_mapper = $self->type_mapper;
    my $seqio =
      Bio::SeqIO->new(-file=> $file,
                      -format => $fmt);
    while (my $seq = $seqio->next_seq()) {
	$unflattener->unflatten_seq(-seq=>$seq,
				    -use_magic=>1);
	$type_mapper->map_types_to_SO(-seq=>$seq);
	my $outio = Bio::SeqIO->new( -format => 'chaos');
	$outio->write_seq($seq);
	my $stag = $outio->handler->stag;
	$self->init_from_stag($stag);
    }
    $self->name_all_features;
    return;
}

# --- turns object into stag document ---

=head2 stag

 Usage   - my $chaos_node = $cg->stag;
 Returns - L<Data::Stag>
 Args    -

Generates a L<Data::Stag> object conforming to Chaos-XML dtd

=cut

sub stag {
    my $self = shift;
    my $W = Data::Stag->makehandler;
    $self->fire_events($W);
    return $W->stag;
}

sub chaos_flavour {
    my $self = shift;
    $self->{_chaos_flavour} = shift if @_;
    return $self->{_chaos_flavour} || 'chaos';
}

sub metadata {
    my $self = shift;
    $self->{_metadata} = shift if @_;
    return $self->{_metadata};
}


sub fire_events {
    my $self = shift;
    my $W = shift;

    my $t = time;
    my $ppt = localtime($t);
    my $prog = $0;
    chomp $prog;

    my @meta = $self->metadata ? ($self->metadata) : ();
    $W->start_event('chaos');
    $W->event(chaos_metadata=>[
			       [chaos_version=>1],
			       [chaos_flavour=>$self->chaos_flavour],
			       @meta,
			       
			       [feature_unique_key=>'feature_id'],
			       [equiv_chado_release=>'chado_1_01'],
			       
			       [export_unixtime=>$t],
			       [export_localtime=>$ppt],
			       [export_host=>$ENV{HOST}],
			       [export_user=>$ENV{USER}],
			       [export_perl5lib=>$ENV{PERL5LIB}],
			       [export_program=>$prog],
			      ]
	     );
    my $g = $self->graph;

    # unordered; features followed by frs
    my $done_idx = {};
    my @ufeats = @{$self->unlocalised_features};
    $self->fire_feature_event($W, $_, $done_idx) foreach @ufeats;
    
    $W->end_event('chaos');
}

sub fire_feature_event {
    my $self = shift;
    my $W = shift           || $self->freak("no writer");
    my $f = shift           || $self->freak("no feature"); 
    my $done_idx = shift    || $self->freak("no index of done features");
    my $fid = $f->get_feature_id;
    return if $done_idx->{$fid};

    my $g = $self->graph;
    my @in_edges = $g->edges_from($fid); # object FRs
    my @frs = ();
    while (my $edge = shift @in_edges) { 
        my ($subject_id,$object_id) = @$edge;
	my $type = $g->get_edge_attribute(@$edge,'type');
	my $rank = $g->get_edge_attribute(@$edge,'rank');

        if (!$self->feature_idx->{$object_id}) {
            $f->add_featureprop([[type=>'comment'],[value=>"this feature has a parent in another subgraph; there will be a trailing object_id=$object_id"]]);
            # this is the case for AceView worm models and
            # dicistronic genes where exons are shared across genes
        }
        else {
            # objects must be written before subjects
            $self->fire_feature_event($W,
                                      $self->feature_idx->{$object_id},
                                      $done_idx);
        }
	# no point carrying on, redundant tree traversal
	return if $done_idx->{$fid};
	push(@frs,
	     [feature_relationship=>[
				     [subject_id=>$subject_id],
				     [object_id=>$object_id],
				     [type=>$type],
				     [rank=>$rank],
				    ]]);
    }
    return if $done_idx->{$fid};

    $W->event(feature=>$f->data);
    $W->event(@$_) foreach @frs;
    $done_idx->{$fid} = 1;
    my @nextfs = 
      (@{$self->get_features_on($f)},
       @{$self->get_features_contained_by($f)});
#    print "$fid has the following: @nextfs\n";
    $self->fire_feature_event($W, $_, $done_idx) foreach @nextfs;
    return;
}


sub init_mldbm {
    my $self = shift;
    require "MLDBM.pm";
    import("MLDBM", qw(DB_File Storable));
    return;
}


# deprecated?
sub next_idn {
    my $self = shift;
    $self->{_next_idn} = shift if @_;
    return $self->{_next_idn};
}


# deprecated?
sub generate_new_feature_id {
    my $self = shift;
    my $prefix = shift || 'feature';
    my $feature_id;
    my $idn = $self->{_next_idn} || 0;
    my $fidx = $self->feature_idx;
    while (!$feature_id) {
	$idn++;
	unless ($fidx->{"$prefix-$idn"}) {
	    $feature_id = "$prefix-$idn";
	}
    }
    $self->{_next_idn} = $idn;
    return $feature_id;
}


=head2 unflattener

  Usage   - my $unf = $cg->unflattener;
  Usage   - $cg->unflattener(Bio::SeqFeature::Tools::Unflattener->new);
  Returns - L<Bio::SeqFeature::Tools::Unflattener>
  Args    - L<Bio::SeqFeature::Tools::Unflattener> (OPTIONAL)

gets/sets the object that the CG will use for unflattening genbank
seqs. See bioperl docs for details

=cut

sub unflattener {
    my $self = shift;
    $self->{_unflattener} = shift if @_;
    if (!$self->{_unflattener} ) {
	$self->load_module("Bio::SeqFeature::Tools::Unflattener");
	$self->{_unflattener} =
	  Bio::SeqFeature::Tools::Unflattener->new;
    }
    return $self->{_unflattener};
}

=head2 type_mapper

  Usage   - my $unf = $cg->type_mapper;
  Usage   - $cg->type_mapper(Bio::SeqFeature::Tools::Type_mapper->new);
  Returns - L<Bio::SeqFeature::Tools::Type_mapper>
  Args    - L<Bio::SeqFeature::Tools::Type_mapper> (OPTIONAL)

gets/sets the object that the CG will use for mapping genbank types to
SO. See bioperl docs for details

=cut

sub type_mapper {
    my $self = shift;
    $self->{_type_mapper} = shift if @_;
    if (!$self->{_type_mapper} ) {
	$self->load_module("Bio::SeqFeature::Tools::TypeMapper");
	$self->{_type_mapper} =
	  Bio::SeqFeature::Tools::TypeMapper->new;
    }
    return $self->{_type_mapper};
}



=head2 feature_idx

 Usage   - my $f = $cg->feature_idx->{$feature_id}
 Returns - hashref, keyed by feature ID
 Args    - 

index hash for looking up feature stag nodes by ID

The ID is the value of feature/feature_id in the chaos-xml structure

the hash has values that are chaos L<Data::Stag> nodes, and can be
accessed using normal stag methods/functions

=cut

sub feature_idx {
    my $self = shift;
    $self->{_feature_idx} = shift if @_;
    return $self->{_feature_idx};
}


=head2 get_feature

 Usage   - my $f = $cg->get_feature($feature_id)
 Returns - L<Data::Stag> conforming to Chaos-xml feature element
 Args    - id string

look up feature stag nodes by ID

The ID is the value of feature/feature_id in the chaos-xml structure

The returned L<Data::Stag> node is a feature node/element, and can be
accessed using normal stag methods/functions

=cut

sub get_feature {
    my $self = shift;
    my $fid = shift;
    return $self->{_feature_idx}->{$fid};
}

# relationship graph
# equiv to chaos/chado feature_relationship graph

=head2 graph

  Usage   - my $graph = $cg->graph;
  Returns - L<Graph>
  Args    -

gets/sets the L<Graph> object which corresponds to the graph defined
by the chado/chaos feature_relationships. The sink of the graph is the
root (eg genes), the source of the graph is the leaves (eg exons).

The graph labels contain the relationship type and the rank

You should not need to manipulate the Graph object directly - this is
mostly used internally, but is made public to allow inspection of the
graph using the native L<Graph> methods

=cut

sub graph {
    my $self = shift;
    $self->{_graph} = shift if @_;
    return $self->{_graph};
}

=head2 graph

  Usage   - my $locgraph = $cg->locgraph;
  Returns - L<Graph>
  Args    -

gets/sets the L<Graph> object which corresponds to the graph defined
by the chado/chaos featurelocs. 

The graph labels contain the nbeg, nend, rank, group

You should not need to manipulate the L<Graph> object directly - this is
mostly used internally, but is made public to allow inspection of the
graph using the native L<Graph> methods

=cut

sub locgraph {
    my $self = shift;
    $self->{_locgraph} = shift if @_;
    return $self->{_locgraph};
}


=head2 add_feature

 Usage   - $cg->add_feature($f);
 Returns -
 Args    - L<Data::Stag> [feature node]

=cut

sub add_feature {
    my $self = shift;
    my $feature = shift;
    my $fid = $feature->get_feature_id;

    my $verbose = $self->verbose;

    $self->graph->add_vertex($fid);
    $self->feature_idx->{$fid} = $feature;
    my @flocs = $feature->get_featureloc;
    foreach my $floc (@flocs) {
	$self->add_featureloc($fid, $floc);
    }
    $self->debug("added feature: $fid") if $verbose;
    return 1;
}

=head2 add_featureloc

 Usage   - $cg->add_featureloc($fl);
 Returns -
 Args    - L<Data::Stag> [featureloc node]

=cut

sub add_featureloc {
    my $self = shift;
    my $fid = shift;
    my $floc = shift;
    my $lg = $self->locgraph;

    my $verbose = $self->verbose;

    my $src_fid = $floc->get_srcfeature_id;
    $lg->add_edge($fid,$src_fid);
    $self->debug("added featureloc: $fid to $src_fid") if $verbose;
    return;
}


=head2 replace_featureloc

  Usage   - $cg->replace_featureloc($f,$fl_old,$fl_new);
  Returns -
  Args    - feature L<Data::Stag> [feature node]
            old loc L<Data::Stag> [featureloc node]
            new loc L<Data::Stag> [featureloc node]

=cut

sub replace_featureloc {
    my $self = shift;
    my $feature = shift;
    my $old_floc = shift;
    my $new_floc = shift;

    my $verbose = $self->verbose;

    my $lg = $self->locgraph;
    my $fid = $feature->get_feature_id;

    my $old_src_fid = $old_floc->get_srcfeature_id;
    my $new_src_fid = $new_floc->get_srcfeature_id;

    my @new_e = ($fid,$new_src_fid);
    $lg->delete_edge($fid, $old_src_fid);
    $lg->add_edge(@new_e);
    foreach ($new_floc->kids) {
        next unless $_->isterminal;
        next if $_->name eq 'srcfeature_id';
        $lg->set_edge_attribute(@new_e,
                                $_->name,
                                $_->data);
    }
    $feature->set_featureloc($new_floc->data);
    $self->debug("replaced featureloc: $fid to $new_src_fid") if $verbose;
    return;
}

=head2 add_feature_relationship

 Usage   - $cg->add_feature_relationship($f);
 Returns -
 Args    - L<Data::Stag> [feature_relationship node]

=cut

sub add_feature_relationship {
    my $self = shift;
    my $fr = shift;
    my $g = $self->graph;
    my %frh = $fr->pairs;
    # FROM subject TO object
    # (this means subject is source, object is sink)
    my @edge = ($frh{subject_id},
                $frh{object_id});

    my $verbose = $self->verbose;

    if (!$edge[0] || !$edge[1]) {
	$self->freak("bad feature_rel", $fr);
    }
    $g->add_edge(@edge);
    $g->set_edge_attribute(@edge,
                           "type",
                           $frh{type} || '');
    $g->set_edge_attribute(@edge,
                           "rank",
                           $frh{rank} || '0');

    $self->debug("added fr: @edge") if $verbose;
    return 1;
}



=head2 top_features

  Usage   - my $features = $cg->top_features;
  Returns - listref of L<Data::Stag> feature nodes
  Args    -

returns features at the root of the feature graph (typically genes,
but also "simple" features that are not attached other features via
feature_relationships, such as SNPs, contigs, etc

Formally, a feature F is a top feature if there is no
feature_relationship R with R.subject_id = F

=cut

sub top_features {
    my $self = shift;
    my $g = $self->graph;
    my $fidx = $self->feature_idx;
    my @fl = ();
    while (my($fid, $f) = each %$fidx) {
	if (!$g->edges_from($fid)) {
	    push(@fl, $f);
	}
    }
    return \@fl;
}


=head2 leaf_features

  Usage   - my $features = $cg->leaf_features;
  Returns - listref of L<Data::Stag> feature nodes
  Args    -

returns features at the leaves of the feature graph (with gene model
subgraphs, these may be exons and polypeptides - or "simple" features
that are not attached other features via feature_relationships, such
as SNPs, contigs, etc)

Formally, a feature F is a leaf feature if there is no
feature_relationship R with R.object_id = F

=cut

sub leaf_features {
    my $self = shift;
    my $g = $self->graph;
    my $fidx = $self->feature_idx;
    my @fl = ();
    while (my($fid, $f) = each %$fidx) {
	if (!$g->edges_to($fid)) {
	    push(@fl, $f);
	}
    }
    return \@fl;
}


=head2 unlocalised_features

  Usage   - my $topfs = $cg->unlocalised_features;
  Synonym - unlocalized_features
  Returns - listref of L<Data::Stag> feature nodes
  Args    -

returns features at the root of the featureloc graph, ie unlocalised
features.

Formally, a feature F is unlocalised if it contains no featurelocs


=cut

sub unlocalised_features {
    my $self = shift;
    my $lg = $self->locgraph;
    my $fidx = $self->feature_idx;
    my @fl = ();
    while (my($fid, $f) = each %$fidx) {
	if (!$lg->edges_from($fid)) {
	    push(@fl, $f);
	}
    }
    return \@fl;
}
*unlocalized_features = \&unlocalised_features;

=head2 top_unlocalised_features

  Usage   - my $topfs = $cg->top_unlocalised_features;
  Synonym - top_unlocalized_features
  Returns - listref of L<Data::Stag> feature nodes
  Args    -

Returns the intersection of the set of all unlocalised features and
all top features

=cut

sub top_unlocalised_features {
    my $self = shift;
    my $g = $self->graph;
    my $lg = $self->locgraph;
    my $fidx = $self->feature_idx;
    my @fl = ();
    while (my($fid, $f) = each %$fidx) {
	if (!$g->edges_from($fid) &&
	    !$lg->edges_from($fid)) {
	    push(@fl, $f);
	}
    }
    return \@fl;
}
*top_unlocalized_features = \&top_unlocalised_features;


=head2 get_features_on

  Usage   - my $features = $cg->get_features_on($contig_feature)
  Returns - listref of L<Data::Stag> feature nodes
  Args    - L<Data::Stag> feature node OR id string

all features DIRECTLY localised to a particular feature

=cut

sub get_features_on {
    my $self = shift;
    my $srcf = shift;

    my $srcfid = ref($srcf) ? $srcf->get_feature_id : $srcf;
    my $lg = $self->locgraph;

    my @edges = $lg->edges_to($srcfid);
    my @located_fids = ();
    while (my $edge = shift @edges) {
	push(@located_fids, $edge->[0]);
    }
    my $fidx = $self->feature_idx;
    return [map {$fidx->{$_}} @located_fids];
}


=head2 get_features_contained_by

  Usage   - my $transcripts = $cg->get_features_contained_by($gene_feature)
  Returns - listref of L<Data::Stag> feature nodes
  Args    - L<Data::Stag> feature node

all features contained by another feature, where containment is
defined by any feature_relationship, with the container being the
object_id and the containee being the
subject_id. feature_relationship.type is ignored

=cut

sub get_features_contained_by {
    my $self = shift;
    my $f = shift;

    my $g = $self->graph;

    my @edges = $g->edges_to($f->get_feature_id);
    my @contained_fids = ();
    while (my $edge = shift @edges) {
	push(@contained_fids, $edge->[0]);
    }
    my $fidx = $self->feature_idx;
    return [map {$fidx->{$_}} @contained_fids];
}

=head2 get_features_containing

  Usage   - my $transcripts = $cg->get_features_containing($exon_feature)
  Returns - listref of L<Data::Stag> feature nodes
  Args    - L<Data::Stag> feature node

all features containing by another feature, where containment is
defined by any feature_relationship, with the container being the
object_id and the containee being the
subject_id. feature_relationship.type is ignored

=cut

sub get_features_containing {
    my $self = shift;
    my $f = shift;

    my $g = $self->graph;

    my @edges = $g->edges_from($f->get_feature_id);
    my @container_fids = ();
    while (my $edge = shift @edges) {
	push(@container_fids, $edge->[1]);
    }
    my $fidx = $self->feature_idx;
    return [map {$fidx->{$_}} @container_fids];
}


=head2 get_all_contained_features

  Usage   - my $features = $cg->get_all_contained_features($gene_feature)
  Returns - listref of L<Data::Stag> feature nodes
  Args    - L<Data::Stag> feature node

As get_features_contained_by, but performs the transitive closure - eg
for a gene will fetch transcriipts, and the transcripts subfeatures
(exons) into one flat list

=cut

sub get_all_contained_features {
    my $self = shift;
    my $top = shift || $self->freak("requires parameter: top [feature]");
    my $topfid = $top->get_feature_id;
    my $fidx = $self->feature_idx;

    my $verbose = $self->verbose;

    $self->debug("get_all_contained_features $topfid") if $verbose;

    my $iterator = $self->feature_iterator($topfid);
    my @cfids = ();
    my %got_idh = ();
    while (my $fid = $iterator->next_vertex) {
        $self->debug("iterator; next=$fid") if $verbose;
        next if $got_idh{$fid};
        $got_idh{$fid} = 1;
	push(@cfids, $fid) unless $fid eq $topfid;
    }
    return [map {$fidx->{$_}} @cfids];
}



=head2 feature_relationships_for_subject

 Usage   - $frs = $cg->get_feature_relationships_for_subject($exon_id);
 Returns - listref of L<Data::Stag> feature_relationship nodes
 Args    - L<Data::Stag> feature node OR id string

find the feature_relationship nodes with a particular subject_id

=cut

sub feature_relationships_for_subject {
    my $self = shift;
    my $f = shift;
    my $fid = ref($f) ? $f->get_feature_id : $f;

    my $g = $self->graph;

    my @edges = $g->edges_from($fid);
    my @frs = ();
    while (my $edge = shift @edges) {
	if (!$edge->[0] || !$edge->[1]) {
	    $self->freak("bad edge: [@$edge]", $f);
	}
	my $type = $g->get_edge_attribute(@$edge,'type');
	push(@frs, 
	     Data::Stag->new(feature_relationship=>[
						    [subject_id=>$edge->[0]],
						    [object_id=>$edge->[1]],
						    [type=>$type],
						   ]));
    }
    return [@frs];
}


sub get_floc {
    my $self = shift;
    my $f = shift;

    my @flocs = $f->get_featureloc;
    if (@flocs > 1) {
	@flocs = grep {!$_->get_rank && !$_->get_locgroup} @flocs;
	if (@flocs > 1) {
	    $self->freak("invalid flocs", @flocs);
	}
    }
    return shift @flocs;
}


=head2 make_gene_islands

 Usage   - my $contigs = $cg->make_gene_islands;
 Returns - listref of L<Data::Stag> feature nodes
 Args    -

create a contig feature for every gene, and transform the gene and all
the subfeatures of gene onto that contig

=cut

sub make_gene_islands {
    my $self = shift;
    my @args = @_;
    my $fs = $self->top_features;
    my @islands = ();
    foreach my $f (@$fs) {
	my $type = $f->get_type;
	$self->freak("no type", $f) unless $type;
	next unless $f->get_type eq 'gene';
	my $island = $self->make_island($f, @args);
	push(@islands, $island);
    }
    return \@islands;
}

# generates an island contig around a feature $f and transforms the
# coordinates to the contig
sub make_island {
    my $self = shift;
    my $f = shift;

    my $verbose = $self->verbose;
    my ($left, $right) = @_;
    if (!$left) {
	$left = 0;
    }
    if (!$right) {
	$right = $left;
    }
    my $floc = $self->get_floc($f);
    if (!$floc) {
	$self->freak("No featureloc", $f);
    }
    my $src_fid = $floc->get_srcfeature_id;
    my $srcf = $self->get_feature($src_fid);
    my $strand = $floc->get_strand;

    my $nbeg = $floc->get_nbeg - $left * $strand;
    my $nend = $floc->get_nend + $right * $strand;
#    my $island_id = $self->generate_new_feature_id('contig');
    my $island_id = "contig:$src_fid:$nbeg:$nend";

    my $island_name = 'contig-'.$f->get_name.'-'.$left.'-'.$right;
    my $island_uniquename = 'contig-'.$f->get_uniquename.'-'.$left.'-'.$right; 
    $self->debug("making island $island_name") if $verbose;

    my $island =
      $self->new_feature(
			 feature_id=>$island_id,
			 name=>$island_name,
			 uniquename=>$island_uniquename,
			 type=>'contig',
			 featureloc=>[
				      nbeg=>$nbeg,
				      nend=>$nend,
                                      strand=>$strand,
				      srcfeature_id=>$src_fid,
				     ],
			);
    $self->debug("deriving residues $island_name") if $verbose;
    $self->derive_residues($island);
    $self->add_feature($island);
    $self->debug("performing main loctransform $island_name") if $verbose;
    $self->loctransform($f,
			$island);
    my $children = $self->get_all_contained_features($f);
    # replicate feature and add to subhraph
    # (we wish to replicate because a feature can be
    #  shared between graphs and we want to do loctransforms
    #  on the features on a per-subgraph basis)
    $children = [map {$_->duplicate} @$children];

    $self->debug("performing subfeatures loctransform $island_name [total %s]",
          scalar(@$children)) if $verbose;
    foreach my $child (@$children) {
	$self->loctransform($child, $island);
    }
    $self->debug("creating new chaos graph for island $island_name") if $verbose;
    my $C = $self->new; # create a new subgraph
    my @feats = ($srcf, $island, $f, @$children);
    foreach my $subf (@feats) {
	$C->add_feature($subf);
	my $frs = $self->feature_relationships_for_subject($subf);
	$C->add_feature_relationship($_) foreach @$frs;
    }
    $self->debug("new chaos graph created for island $island_name") if $verbose;
    return $C;
}


=head2 derive_residues

 Usage   - my $ok = $cg->derive_residues($feature);
 Returns - sequence string
 Returns - L<Data::Stag> feature node

splices out the residues from the srcfeature and sets the
feature/residues element (does not return the actual residues)

=cut

sub derive_residues {
    my $self = shift;
    my $feature = shift;
    my $res;
    if ($self->is_spliced($feature)) {
        $self->freak('not yet');
    }
    else {
        my @flocs = $feature->get_featureloc;
	if (!@flocs) {
	    $self->freak("feature is not located, can't derive residues",
			 $feature);
	}
        @flocs = grep {!$_->get_rank} @flocs;
        $self->freak unless @flocs;
        my @resl =
          map {
	      my $src_fid = $_->get_srcfeature_id;
              my $srcf = $self->get_feature($src_fid);
	      if (!$srcf) {
		  $self->freak("no source feature for $src_fid in feature",
			       $feature);
	      }
              my $srcres = $srcf->get_residues;
	      if (!$srcres) {
		  $self->freak("feature $src_fid has no residues", $srcf);
	      }
              $self->cutseq($srcres, $_->get_nbeg, $_->get_nend);
          } @flocs;
        $res = shift @resl;
        if (@resl) {
            foreach (@resl) {
                if ($_ ne $res) {
                    $self->freak("$_ ne $res");
                }
            }
        }
        
    }
    $self->freak("cannot derive residues", $feature) unless defined $res;
    $feature->set_residues($res);
    return 1;
}

sub cutseq {
    my $self = shift;
    my $res = shift;
    my $nbeg = shift;
    my $nend = shift;
    if ($nbeg <= $nend) {
        return substr($res, $nbeg, $nend-$nbeg);
    }
    else {
        my $cut = substr($res, $nend, $nbeg-$nend);
        $cut = $self->revcomp($cut);
        return $cut;
    }
}

sub revcomp {
    my $self = shift;
    my $res = shift;
    $res =~ tr/acgtrymkswhbvdnxACGTRYMKSWHBVDNX/tgcayrkmswdvbhnxTGCAYRKMSWDVBHNX/;
    return scalar(CORE::reverse($res));
}


=head2 loctransform

 Usage   - $cg->loctransform($gene,$contig);
 Args    - source L<Data::Stag> feature node
           target L<Data::Stag> feature node

replaces the featureloc(s) of a feature with new featureloc(s)
relative to the target feature - eg going from chromosome to contig

=cut

sub loctransform {
    my $self = shift;
    my $sfeature = shift;                    # source  (eg gene)
    my $tfeature = shift;                    # target  (eg contig)

    my $verbose = $self->verbose;

    # get source and target feature locations;
    # any feature can have >1 flocs (differentiated by rank, group)
    # (usually there will be just 1 each)
    my @sflocs = $sfeature->get_featureloc;
    my @tflocs = $tfeature->get_featureloc;

    $self->debug("  loctransform srcs:%s targets:%s",
                 scalar(@sflocs),scalar(@tflocs)) if $verbose;
    # the source and target locations we use to actually transform
    my $sfloc;
    my $tfloc;
    my $ssrc_fid;
    my $tsrc_fid;

    my $already_transformed;
    # ASSERTION:
    # forall (@sflocs, @tflocs)
    #     there exists exactly one pair ($sfloc, $tfloc)
    #     such that $sfloc and $tfloc share the same srcfeature_id
    #
    # this pair is the source and target flocs that will be used in
    # the location transform
    foreach my $sflocI (@sflocs) {
        my $ssrc_fidI = $sflocI->get_srcfeature_id;
        if ($ssrc_fidI eq $tfeature->get_feature_id) {
            $already_transformed = 1;
            last;
        }
        foreach my $tflocI (@tflocs) {
            my $tsrc_fidI = $tflocI->get_srcfeature_id;
            # intersection
            if ($ssrc_fidI eq $tsrc_fidI) {
                if ($sfloc || $tfloc) {
                    $self->freak("CONFLICT: >1 pair [$ssrc_fid, $tsrc_fid]",
                                 $sfloc, $tfloc);
                }
                $sfloc = $sflocI;
                $tfloc = $tflocI;
                $ssrc_fid = $ssrc_fidI;
                $tsrc_fid = $tsrc_fidI;
            }
        }
    }

    if ($already_transformed) {
        # nothing to be done
        return;
    }

    # ASSERTION (see above) - at least 1 pair
    if (!($sfloc || $tfloc)) {
        $self->freak("NO LOC PAIR FOUND",
                     @sflocs, @tflocs,$sfeature,$tfeature);
    }


    # s: source
    # t: target

    my $snbeg = $sfloc->get_nbeg;
    my $snend = $sfloc->get_nend;
    my $srank = $sfloc->get_rank;
    my $sstrand = $sfloc->get_strand;
    my $tnbeg = $tfloc->get_nbeg;
    my $tnend = $tfloc->get_nend;
    my $tstrand = $tfloc->get_strand;

    my $tfid = $tfeature->get_feature_id;
    if (!$tfid) {
	$self->freak("NO FEATURE_ID", $tfeature);
    }

    $snbeg = ($snbeg - $tnbeg) * $tstrand;
    $snend = ($snend - $tnbeg) * $tstrand;

    $self->debug("    new floc: %s..%s on %s",
                 $snbeg,$snend,$tfid)
      if $verbose;

    my $nu_sfloc =
      $self->new_featureloc(srcfeature_id=>$tfid,
                            nbeg=>$snbeg,
                            nend=>$snend,
                            strand=>$sstrand,
                            rank=>$srank,
                           );
    $self->replace_featureloc($sfeature, $sfloc, $nu_sfloc);
    $self->debug("  performed loctransform srcs:%s targets:%s",
                scalar(@sflocs),scalar(@tflocs)) if $verbose;
    return;
}

sub history_log {

}

sub new_feature {
    my $self = shift;
    return
      Data::Stag->unflatten(feature=>[@_]);
}

sub new_featureloc {
    my $self = shift;
    return
      Data::Stag->unflatten(featureloc=>[@_]);
}

our %SPLICEDF =
  (mRNA=>1);
sub is_spliced {
    my $self = shift;
    my $feature = shift;
    my $type = $feature->get_type;
    return $SPLICEDF{$type} || 0;
}

sub iterate {
    my $self = shift;
    my $G = shift;
    my $v = shift;
    my $func = shift;
    my $iterator = $self->iterator($G, $v);
    while (my $next_v = $iterator->next_vertex) {
	$func->($next_v);
    }
}

sub iterator {
    my $self = shift;
    return Iterator->new(@_);
}

sub feature_iterator {
    my $self = shift;
    return Iterator->new($self->graph, @_);
}


=head2 get_features_by_type

 Usage   - my $exons = $cg->get_features_by_type('exon');
 Returns - listref of L<Data::Stag> feature nodes
 Args    - type string

gets features by type (exact - does not traverse ontology graph)

=cut

sub get_features_by_type {
    my $self = shift;
    my $type = shift;
    my $fidx = $self->feature_idx;
    my @fs = grep {$_->get_type eq $type} values %$fidx;
    return [@fs];
}

=head2 get_features

 Usage   - my $features = $cg->get_features;
 Returns - listref of L<Data::Stag> feature nodes
 Args    - none

returns all features

=cut

sub get_features {
    my $self = shift;
    my $fidx = $self->feature_idx;
    return [values %$fidx];
}

sub validate {
    my $self = shift;
    my $W = shift;
    my $G = $self->graph;
    my $fidx = $self->feature_idx;
    my $vertices = $G->vertices;
    $W->start_event('chaos_validation');
    my @missing_fids = ();
    my @errs = ();
    foreach my $v (@$vertices) {
	if (exists $fidx->{$v}) {
	}
	else {
	    $W->event(missing_feature=>$v);
	    push(@missing_fids, $v);
	}
    }
    if (@missing_fids) {
	push(@errs, "Missing feature_ids: @missing_fids");
    }
    my $features = $self->get_features;
    foreach my $f (@$features) {
	my $name = $f->get_name;
	my $res = $f->get_residues;
	my @flocs = $f->get_featureloc;
	if ($res && scalar(@flocs)) {
	    my $implicit_res = $self->derive_residues($f);
	    if ($res ne $implicit_res) {
		$W->event(residues_conflict=>$name);
		push(@errs, "residues $name");
	    }
	}
    }
    $W->end_event('chaos_validation');
    return @errs;
}


=head2 name_all_features

 Usage   - $cg->name_all_features
 Returns -
 Args    -

makes sure all feature have names. will not affect features that
already have names

sets both feature/name and feature/uniquename

=cut

sub name_all_features {
    my $self = shift;
    my $basename = shift;

    my %global_id_by_type = ();   # for unnamed top features

    my $topfs = $self->top_features;
    foreach my $topf (@$topfs) {
	my $childfs = $self->get_all_contained_features($topf);
	
	my $tname = $topf->get_name;
	if (!$tname) {
	    my $type = $topf->get_type;
	    $global_id_by_type{$type} = 0 unless $global_id_by_type{$type};
	    my $id = ++$global_id_by_type{$type};
	    $tname = "$type$id";
	    if ($basename) {
		$tname = "$basename-$tname";
	    }
	    $topf->set_name($tname);
#	    $topf->set_uniquename($tname);
	}
	my %id_by_type = ();      # unique within a topfeature
    
	foreach my $cf (@$childfs) {
	    my $type = $cf->get_type;
	    $id_by_type{$type} = 0 unless $id_by_type{$type};
	    my $id = ++$id_by_type{$type};
	    my $name = "$tname-$type-$id";
	    $cf->set_name($name);
	    $cf->set_uniquename($name);
	}
    }
    return;
}

sub asciitree {
    my $self = shift;
    my $containers = $self->unlocalised_features;
    my $fidx = $self->feature_idx;
    foreach my $f (@$containers) {
	$self->asciifeature($f, 0);
    }
}

sub asciifeature {
    my $self = shift;
    my $f = shift;
    my $indent = shift || 0;

    my @flocs = $f->get_featureloc;
    printf("%s%s %s \"%s\" %s\n",
	   ' '  x $indent,
	   $f->get_type,
	   $f->get_feature_id,
	   $self->get_feature_shortlabel($f->get_feature_id),
	   join(";",
		map {
		    sprintf("%s->%s on %s",
			    $_->get_nbeg, $_->get_nend,
			    $self->get_feature_shortlabel($_->get_srcfeature_id))
		} @flocs),
	  );
    my $cfeats = $self->get_features_contained_by($f);
    foreach my $subf (@$cfeats) {
	$self->asciifeature($subf, $indent+1);
    }
    my $lfeats = $self->get_features_on($f);
    foreach my $subf (@$lfeats) {
	my $parents = $self->get_features_containing($subf);
	next if @$parents; # roots only
	printf("%s[anchors]\n",
	       ' ' x ($indent+1));
	$self->asciifeature($subf, $indent+2);
    }
}

sub get_feature_shortlabel {
    my $self = shift;
    my $fid = shift;
    my $fidx = $self->feature_idx;
    my $f = $fidx->{$fid};
    return '?' unless $f;
    my $name = $f->get_name;
    return $name if $name;
    return $fid;
}

sub debug {
    my $self= shift;
    return unless $self->verbose;
    my $fmt = shift;
    my $t = time;
    my $lt = localtime $t;
    print STDERR "# "; 
    printf STDERR ($fmt, @_);
    print STDERR " [$t] $lt\n";
}


1;

package Iterator;

sub new {
    my $self = shift;
    my $G = shift;    # graph or array of graphs
    my $v = shift;
    my $dir = shift || 'out';
    unless (ref($G) eq 'ARRAY') {
	$G = [$G];
    }

    # if vertices not supplied by user, choose a random one
    if (!$v) {
	($v) = map {@{$_->vertices}} @$G;
    }
    my $depth = 0;

    # nodes to be explored;
    #  each node is an arrayref [$depth,$subj,$obj,$rank]
    my @nodes = ();
    
    my $closure = sub {
	my $meth = shift;
	if ($meth eq 'next_vertex') {
	    my @edges;
	    my @all_child_nodes = ();
	    foreach my $g (@$G) {
		if ($dir eq 'out') {
                    # sink-to-source (root-to-leaf)
		    @edges = $g->edges_to($v);
		}
		else {
                    # source-to-sink (leaf-to-root) DEFAULT
		    @edges = $g->edges_from($v);
		}
		my @child_nodes = ();
		while (my $edge = shift @edges) {
		    my ($subj,$obj) = @$edge;
		    my $rank = $g->get_edge_attribute($subj,$obj,'rank');
		    die("Assertion error") unless defined $rank;
		    push(@child_nodes, [$depth+1, $subj, $obj, $rank]);
		}
		@child_nodes = sort { $a->[3] <=> $b->[3] } @child_nodes;
                # default: depth-first
		push(@all_child_nodes, @child_nodes);
	    }
	    push(@nodes, @all_child_nodes);
	    my $nextnode = shift @nodes;
	    if (!$nextnode) {
		$depth = -1;
		return;
	    }
	    $depth = $nextnode->[0];
            # out: use subj
            # in:  use obj
	    $v = $dir eq 'out' ? $nextnode->[1] : $nextnode->[2];
	    return $v;
	}
	elsif ($meth eq 'depth') {
	    return $depth;
	}
	else {
	    $self->freak("cannot call method \"$meth\" on an iterator");
	}
    };
    bless $closure, 'Iterator';
    return $closure;
}

sub next_vertex { &{shift @_}('next_vertex')}
sub depth { &{shift @_}('depth')}


1;