package Bio::Graphics::Glyph::pairplot;

# Triangle plot for showing pairwise quantitative relationships.
# Similar to pairwise_plot, which was originally distributed in GBrowse,
# except that the signal intensity is handled using a callback that takes
# the two features to be compared.

sub my_description {
    return <<END;
This is a triangle plot that displays quantitative relationships between point-like features, 
such as LD between SNPs. The quantitative information is generated by the "score" callback,
which receives two feature objects and returns a floating point number between "min_score"
and "max_score". The bgcolor intensity will be scaled from min to max like graded_segments.

NOTE (2 April 2008): This glyph is not yet fully functional.
END
}

sub my_options {
    {
	score => [
	    'CODEREF',
	    undef,
	    'Pass a coderef to a subroutine with the following signature: sub ($$). The',
	    'two arguments are the first and second feature to compare. Return a floating point',
	    'number to indicate the score at the intersection of these two features.'
	],
	min_score => [
	    'float',
	    0.0,
	    'Minimum possible pairwise score.'
	],
	max_score => [
	    'float',
	    1.0,
	    'Maximum possible pairwise score.'
	],
        angle => [
	    'float',
	    45,
	    'Angle between the side of the triangle and the base, in degrees.'
        ],
    }
}



use strict;
use Math::Trig;

use base 'Bio::Graphics::Glyph::generic';

sub maxdepth {
  my $self = shift;
  my $md   = $self->Bio::Graphics::Glyph::maxdepth;
  return $md if defined $md;
  return 1;
}

# return angle in radians
sub angle {
  my $self  = shift;
  my $angle = $self->{angle} ||= $self->option('angle') || 45;
  $self->{angle} = shift if @_;
  deg2rad($angle);
}

sub slope {
  my $self = shift;
  return $self->{slope} if exists $self->{slope};
  return $self->{slope} = tan($self->angle);
}

sub x2y {
  my $self = shift;
  shift() * $self->slope;
}

sub intercept {
  my $self = shift;
  my ($x1,$x2) = @_;
  my $mid = ($x1+$x2)/2;
  my $y   = $self->x2y($mid-$x1);
  return (int($mid+0.5),int($y+0.5));
}

# height calculated from width
sub layout_height {
  my $self = shift;
  return $self->{height} if exists $self->{height};
  return $self->{height} = $self->x2y($self->width)/2;
}

sub min_score {
    my $self = shift;
    my $min  = $self->option('min_score');
    return   0 unless defined $min;
}
sub max_score {
    my $self = shift;
    my $max  = $self->option('max_score');
    return 1.0 unless defined $max;
}

sub calculate_color {
  my $self = shift;
  my ($s,$rgb) = @_;
  return $self->{colors}{$s} if exists $self->{colors}{$s};
  my $max_score = $self->max_score;
  my $min_score = $self->min_score;
  my $scale     = 255/($max_score-$min_score);
  my $value     = ($s-$min_score) * $scale; # will range from 0 to 255
  return $self->{colors}{$s} = 
    $self->panel->translate_color(map { 255 - $value} @$rgb);
}

sub draw {
  my $self = shift;
  my $gd   = shift;
  my ($left,$top,$partno,$total_parts) = @_;
  my $fgcolor = $self->fgcolor;

  my ($red,$green,$blue) = $self->panel->rgb($self->bgcolor);

  my @points = $self->get_points();
  $gd->line($self->left+$left, $top+1,
	    $self->right+$left,$top+1,
	    $fgcolor);

  my $points = $self->option('point');

  my @parts = sort {$a->left<=>$b->left} $self->parts;
  $_->draw_component($gd,$left,$top-10) foreach @parts;

  # assumption: parts are not overlapping
  if ($points) {
    @points = map { int (($parts[$_]->right+$parts[$_+1]->left)/2)} (0..$#parts-1);
    unshift @points,int($parts[0]->left);
    push @points,int($parts[-1]->right);
  }

  for (my $ia=0;$ia<@parts-1;$ia++) {
    for (my $ib=$ia+1;$ib<@parts;$ib++) {
      my ($l1,$r1,$l2,$r2);
      if (@points) {
	($l1,$r1) = ($points[$ia]+1,$points[$ia+1]-1);
	($l2,$r2) = ($points[$ib]+1,$points[$ib+1]-1);
      } else {
	($l1,$r1) = ($parts[$ia]->left,$parts[$ia]->right);
	($l2,$r2) = ($parts[$ib]->left,$parts[$ib]->right);
      }

      my $intensity = eval{$self->feature->pair_score($parts[$ia],$parts[$ib])};
      warn $@ if $@;
      $intensity    = 1.0 unless defined $intensity;
      my $c         = $self->calculate_color($intensity,[$red,$green,$blue]);

      # left corner
      my ($lcx,$lcy) = $self->intercept($l1,$l2);
      my ($tcx,$tcy) = $self->intercept($r1,$l2);
      my ($rcx,$rcy) = $self->intercept($r1,$r2);
      my ($bcx,$bcy) = $self->intercept($l1,$r2);
      my $poly = GD::Polygon->new();
      $poly->addPt($lcx+$left,$lcy+$top);
      $poly->addPt($tcx+$left,$tcy+$top);
      $poly->addPt($rcx+$left,$rcy+$top);
      $poly->addPt($bcx+$left,$bcy+$top);
      $gd->filledPolygon($poly,$c);
    }
  }
}

sub get_points {
  my $self = shift;
  my @points;
  my @parts = $self->parts;
  return unless @parts;

  for my $g (@parts) {
    push @points,$g->left;
    push @points,$g->right;
  }
  @points;
}

# never allow our internal parts to bump;
sub bump { 0 }

1;

__END__

=head1 NAME

Bio::Graphics::Glyph::pairplot - The "pairwise plot" glyph

=head1 SYNOPSIS

 use Bio::Graphics;

 # create the panel, etc.  See Bio::Graphics::Panel
 # for the synopsis

 # Create one big feature using the PairFeature
 # glyph (see end of synopsis for an implementation)
 my $block = PairFeature->new(-start=>  2001,
 			      -end  => 10000);

 # It will contain a series of subfeatures.
 my $start = 2001;
 while ($start < 10000) {
   my $end = $start+120;
   $block->add_SeqFeature($bsg->new(-start=>$start,
				    -end  =>$end
				   ),'EXPAND');
   $start += 200;
 }

 $panel->add_track($block,
 		   -glyph => 'pairplot',
		   -angle => 45,
		   -bgcolor => 'red',
		   -point => 1,
		  );

 print $panel->png;

 package PairFeature;
 use base 'Bio::SeqFeature::Generic';

 sub pair_score {
   my $self = shift;
   my ($sf1,$sf2) = @_;
   # simple distance function
   my $dist  = $sf2->end    - $sf1->start;
   my $total = $self->end   - $self->start;
   return sprintf('%2.2f',1-$dist/$total);
 }

=head1 DESCRIPTION

This glyph draws a "triangle plot" similar to the ones used to show
linkage disequilibrium between a series of genetic markers.  It is
basically a dotplot drawn at a 45 degree angle, with each
diamond-shaped region colored with an intensity proportional to an
arbitrary scoring value relating one feature to another (typically a
D' value in LD studies).

This glyph requires more preparation than other glyphs.  First, you
must create a subclass of Bio::SeqFeature::Generic (or
Bio::Graphics::Feature, if you prefer) that has a pair_score() method.
The pair_score() method will take two features and return a numeric
value between 0.0 and 1.0, where higher values mean more intense.

You should then create a feature of this new type and use
add_SeqFeature() to add to it all the genomic features that you wish
to compare.

Then add this feature to a track using the pairplot glyph.  When
the glyph renders the feature, it will interrogate the pair_score()
method for each pair of subfeatures.

=head2 OPTIONS

In addition to the common options, the following glyph-specific
options are recognized:

  Option      Description                  Default
  ------      -----------                  -------

  -point      If true, the plot will be         0
              drawn relative to the
              midpoint between each adjacent
              subfeature.  This is appropriate
              for point-like subfeatures, such
              as SNPs.

  -angle      Angle to draw the plot.  Values   45
              between 1 degree and 89 degrees
              are valid.  Higher angles give
              a more vertical plot.

  -bgcolor    The color of the plot.            cyan

=head1 BUGS

Please report them.

=head1 SEE ALSO

L<Bio::Graphics::Panel>,
L<Bio::Graphics::Glyph>,
L<Bio::Graphics::Glyph::arrow>,
L<Bio::Graphics::Glyph::cds>,
L<Bio::Graphics::Glyph::crossbox>,
L<Bio::Graphics::Glyph::diamond>,
L<Bio::Graphics::Glyph::dna>,
L<Bio::Graphics::Glyph::dot>,
L<Bio::Graphics::Glyph::ellipse>,
L<Bio::Graphics::Glyph::extending_arrow>,
L<Bio::Graphics::Glyph::generic>,
L<Bio::Graphics::Glyph::graded_segments>,
L<Bio::Graphics::Glyph::heterogeneous_segments>,
L<Bio::Graphics::Glyph::line>,
L<Bio::Graphics::Glyph::pinsertion>,
L<Bio::Graphics::Glyph::primers>,
L<Bio::Graphics::Glyph::rndrect>,
L<Bio::Graphics::Glyph::segments>,
L<Bio::Graphics::Glyph::ruler_arrow>,
L<Bio::Graphics::Glyph::toomany>,
L<Bio::Graphics::Glyph::transcript>,
L<Bio::Graphics::Glyph::transcript2>,
L<Bio::Graphics::Glyph::translation>,
L<Bio::Graphics::Glyph::triangle>,
L<Bio::Graphics::Glyph::xyplot>,
L<Bio::DB::GFF>,
L<Bio::SeqI>,
L<Bio::SeqFeatureI>,
L<Bio::Das>,
L<GD>

=head1 AUTHOR

Lincoln Stein E<lt>lstein@cshl.edu<gt>.

Copyright (c) 2004 Cold Spring Harbor Laboratory

This package and its accompanying libraries is free software; you can
redistribute it and/or modify it under the terms of the GPL (either
version 1, or at your option, any later version) or the Artistic
License 2.0.  Refer to LICENSE for the full license text. In addition,
please see DISCLAIMER.txt for disclaimers of warranty.

=cut