#!/usr/bin/perl
# -*-Perl-*-
# Last changed Time-stamp: <2007-04-26 14:27:02 xtof>
# $Id$

use Getopt::Long;
use Pod::Usage;
use Data::Dumper;
use strict;
use warnings;
use vars qw/$T0 $T8 $TX $P0 $TINC @FILES $TREEKIN %ENV/;

# defaults for global(s)
$TREEKIN = "$ENV{HOME}/C/treekin/treekin";
$T0 = 0.1;
$T8 = 10;
$TX = -1;
$TINC = 1.02;
$P0 = "";

# process command-line
pod2usage(-verbose => 0)
    unless GetOptions("t0=f"  => \$T0,
		      "t8=f"  => \$T8,
		      "tX=f"  => \$TX,
		      "inc=f" => \$TINC,

		      "p0=s"  => \&process_opt_p0,
		      "man"   => sub{pod2usage(-verbose => 2)},
		      "help"  => sub{pod2usage(-verbose => 1)});

my ($rs, $cs)  = parse_bar_map();
my $p_zeros    = $P0;
my $start_time = $T0;
my $stop_time  = calculate_stop_time($start_time);
my $global_time = 0.0;

# main loop
my $flag = 0;
for (my $file = 0; $file<=$#FILES; $file++) {
  $stop_time = $TX, $flag = 1 if ($#FILES == $file) && ($TX > 0);
  my $command = build_command($p_zeros,
			      $stop_time,
			      $FILES[$file]);
  my ($stop, $densities, $tc) = do_simulation($command, $flag);
  if ($stop == -1) {
    print STDERR "Treekin ERROR at inputfile $FILES[$file]\n";
    last;
  }
  # dump time course to STDOUT
  dump_time_course($stop, $global_time, $tc, $flag);
  $global_time += $stop;
  
  last if $file == $#FILES;
  $p_zeros = remap_densities($densities, $cs->[$file], $cs->[$file+1]);
}

#---
sub dump_time_course {
  my ($stop, $global_time, $tc, $flag) = @_;
  
  # print time course
  foreach my $line (@$tc) {
    if ($line =~ m/^\#/) {
      print $line;
      next;
    }
    
    $line =~ s/^\s+//;
    my @values = split /\s+/, $line;
    $values[0] += $global_time;
    print join(" ", @values), "\n";
  }
  print "&\n";

  return if $flag;
  $stop += $global_time;
  print
      "\@with line\n",
      "\@line on\n",
      "\@line loctype world\n",
      "\@line g0\n",
      "\@line $stop, 1, $stop, 0\n",
      "\@line linewidth .1\n",
      "\@line linestyle 1\n",
      "\@line color 7\n",
      "\@line arrow 0\n",
      "\@line arrow type 0\n",
      "\@line arrow length 1.000000\n",
      "\@line arrow layout 1.000000, 1.000000\n",
      "\@line def\n";
}

#---
sub calculate_stop_time {
  my $start = shift;
  my $time;
  for ($time=$start; $time <= $start+$T8; $time *= $TINC) {;}
  return $time;
}

#---
sub remap_densities {
  my ($densities, $current, $next) = @_;
  my $tuples = [ grep {!     (($_->[0] == -1)
			   || ($_->[1] == -1))} @{zip($current, $next)}];
  # there can be duplicate entries in the barmap file
  $tuples = uniq_tuples($tuples);

  my @p = (0) x ($#$densities);

  for(my $i=0; $i<= $#$tuples; $i++) {
    next if $tuples->[$i]->[0] < 0;
    $p[$tuples->[$i]->[1]] += $densities->[$tuples->[$i]->[0]];
  }

  my $yy = zip(number_sequence(), \@p);
  my $p0 = "";
  my $p0_sum = 0.0;
  for my $tuple (@$yy) {
    next if $tuple->[1] == 0;
    $p0_sum += $tuple->[1];
    $p0 .= " --p0 $tuple->[0]=$tuple->[1]" 
  }

  if ($p0_sum > 1.05) {
    print STDERR "$p0\np0_sum = $p0_sum\n";    
    print STDERR Dumper([grep {!  (($_->[0] == -1)
			        || ($_->[1] == -1))} @$tuples]);
    print STDERR Dumper($yy);
  }

  return $p0;
}

#---
sub uniq_tuples {
  my $tuples = shift;
  my @uniq = ();
  my $h = {map { $_->[0] => $_->[1] } @$tuples};
  foreach my $key (sort {$a <=> $b} keys %$h) {
    push @uniq, [$key, $h->{$key}];
  }
  return \@uniq;
}

#---
# returnes true if first and second of tuple is greater -1
#sub is_defined {
#  return ! (($_[0]->[0] == -1) || ($_[0]->[1] == -1));
#}

#---
sub build_command {
  my ($p_zeros, $stop, $file) = @_;
  my $command = "$TREEKIN -m I";
  $command .= " --tinc=$TINC";
  $command .= " --t0=0.0";
  $command .= " --t8=$stop";
  $command .= $p_zeros;
  $command .= " $file";

  return $command;
}

#---
sub do_simulation {
  my ($command, $flag) = @_;
  print STDERR "$command\n\n";
  my @output = `$command`;

  # capture treekin error
  if ($#output == $[-1) {
    return (-1, undef, undef); 
  }

  # get last line of time course
  $output[-2] = '#' . $output[-2] if $flag == 1;
  my $lastline = $output[-2];
  
  $lastline =~ s/^\s+//;
  # extract values from last line of time course
  my @values = split /\s+/, $lastline;
  # get stop time
  my $stop = $values[0];

  return ($stop, \@values, \@output);
}

#---
sub filter { [map { $_[0]->($_) ? $_ : ()} @{$_[1]}] }

#---
#sub is_defined { return 1 if defined $_->[1]; return 0 }

#---
sub process_opt_p0 {
  my $string = $_[1];
  my ($lmin, $value) = split(/=/, $string);
  $P0 .= " --p0 $lmin=$value";
}

#---
sub parse_bar_map {
  my (@rows, $columns);
  while (<>) {
    s/^\#//, @FILES = split, next if $. == 1;
    chomp;
    
    # smash line into substrings of length 3 and capture only the odd ones
    # trim whitespace and convert empty strings to undef values
    my $row = [map {s/ //g; $_ eq "" ? -1 : $_ } (m/(...)(?:...)?/g)];
    # memorize row in a lol
    push @rows , $row;
    $columns = to_columns($#$row) if $. == 2;
    # memorize columns in a lol
    $columns->(zip(number_sequence(), $row));
  }

  # return lol's of rows and columns
  return (\@rows, $columns->());
}

#---
# zip two lists togather into a list of tuples
# [a1, a2, ...], [b1, b2, ...] -> [[a1, b1], [a2, b2], ...]
sub zip {
  my ($A, $B) = @_;
  $A = array_to_iterator($A) if ref($A) eq 'ARRAY';
  $B = array_to_iterator($B) if ref($B) eq 'ARRAY';
  my @result;
  while (   defined(my $a = $A->())
	 && defined(my $b = $B->())) {
    push @result, [$a, $b];
  }

  return \@result;
}

#---
# unzip a list of tupels into a tuple of lists
# [[a1, b1], [a2, b2], ...] -> [[a1, a2, ..], [b1, b2, ...]]
sub unzip {
  my (@A, @B);
  push(@A, $_->[0]), push(@B, $_->[1]) for @{shift()};

  return [\@A, \@B];
}

#---
# (val1, val2, ...) | [val1, val2, ...] -> {}
# converts an array into an iterator
sub array_to_iterator {
  my @array = (ref $_[0] eq 'ARRAY') ? @{shift()} : @_;

  return sub { shift @array }
}

#---
# n -> { }
# lazy version of infinit list [n..] of integers starting at value n
# implemented as closure
sub number_sequence {
  my $value = defined $_[0] ? $_[0] : 0;

  return sub { $value++ }
}

#---
# incrementally build lol from list of duples, implemented as closure
# [[idx1,val1], [idx2,val2], ...] -> [[values, val1],[values, val2], ...]
sub to_columns {
  my $lol = [ map{[]} $[..$_[0] ]; # empty lol with right dimension
  return sub {
    return $lol unless defined $_[0]; # return result
    map {push @{$lol->[$_->[0]]}, $_->[1]} @{$_[0]}; 
  }
}

=head1 NAME

barmap_simulator - perform dynamic landscape simulation

=head1 SYNOPSIS

barmap_simulator.pl [[-t0 I<FLOAT>] [-t8 I<FLOAT>] [-inc I<FLOAT>]] -p0 I<STRING> barmap-file

=head1 DESCRIPTION

The program writes a shell-script which drives the program C<treekin>
to simulate RNA folding on a dynamic landscape.

=head1 OPTIONS

=over 4

=item B<-t0> I<FLOAT>

Set start time of kinetic simulation to I<FLOAT> (default: 0.1).

=item B<-t8> I<FLOAT>

Set stop time of kinetic simulation to I<FLOAT> (default: 10.0).

=item B<-inc> I<FLOAT>

Set time increment of kinetic simulation to I<FLOAT> (default: 1.02).

=item B<-p0> I<STRING>

Specify initial population for minima on the first energy
landscape. For example if we want to start with 35% of the population
in minimum 3 and 65% in minimum 5 the I<STRING> looks as follows
"3=.35:5=.65" (floats must sum to 1).

=back

=head1 AUTHORS

Christoph Flamm

=head1 BUGS

Please send comments and bug reports to xtof@tbi.univie.ac.at

=cut
    
__END__
# my $a1 = [qw/a b c d e f g/];
# my $b1 = [qw/h i j k l m n/];
# my $c1 = [qw/o p q r s t u/];

# my $a2 = zip_with_natural_numbers($a1, natural_numbers());
# my $b2 = zip_with_natural_numbers($b1, natural_numbers());
# my $c2 = zip_with_natural_numbers($c1, natural_numbers());

# print Dumper($a2);
# print Dumper($b2);

# my $cols = to_cols($#$a2);
# $cols->($a2);
# print Dumper($cols->());
# $cols->($b2);
# $cols->($c2);
# print Dumper($cols->());

#my $a3 = unzip($a2);
#print Dumper($a3);