# creation of a TFBS::MatrixSet object
    # let @list_of_matrix_objects be a list of TFBS::Matrix::* objects

    ###################################
    # Create a TFBS::MatrixSet object:

    my $matrixset = TFBS::MatrixSet->new(); # creates an empty set
    $matrixset->add_Matrix(@list_of_matrix_objects); #add matrix objects to set
    $matrixset->add_Matrix($matrixobj); # adds a single matrix object to set

    # or, same as above:

    my $matrixset = TFBS::MatrixSet->new(@list_of_matrix_objects, $matrixobj);

    ###################################
    #

    my ($self, %args) = @_;
    return TFBS::_Iterator::_MatrixSetIterator->new($self->{matrix_list},
				$args{'-sort_by'},
				$args{'-reverse'}
			       );

    my ($self,%args) = @_; #the rest of @_ goes to _to_seqob;
    my @PWMs;

    # prepare the sequence

    my $seqobj = $self->_to_seqobj(%args);
    $seqobj->{_pdl_matrix} = _seq_to_pdlmatrix($seqobj);

    # prepare the PWMs

    $self->reset();

    while (my $pwm = $self->next() ) {
	push @PWMs,$pwm;
    }

    # do the analysis

    my $hitlist = TFBS::SiteSet->new();

    foreach my $pwm (@PWMs)  {
	my $threshold = ($args{-threshold} or $pwm->{minscore});
	$hitlist->add_siteset($pwm->bsearch(-seqobj=>$seqobj,
					    -threshold =>$threshold ));
    }
    delete $seqobj->{_pdl_matrix};
    return $hitlist;

    my ($self, $cluster, $node) = @_;

    if ($node->is_Leaf()) {
	for (@{$self->{matrix_list}}) {
	    if ($_->ID() eq $node->id()) {
		$cluster->add_matrix($_);
		return;
	    }
	}
    } else {
	$self->_build_cluster($cluster, $_) for ($node->each_Descendent());
    }

    my ($self, $opt, $arg, $list) = @_;

    # Invalid argument
    if (not defined($list)) {
	$error = "Invalid argument: $opt\n";
	return FALSE;     
    }

    # Valid flag or switch.
    return TRUE if (not scalar(@$list)); 


    # Valid switch, check the argument
    for (@$list) {
	return TRUE if ($arg eq $_) ;
    }

    # Valid switch, invalid argument
    $error = "$arg is invalid argument to $opt";
    return FALSE;

    my ($self, %args) = @_;
    my $PWM_SEARCH = '/home/httpd/cgi-bin/CONSITE/bin/pwm_searchPFF';

    # DIRTY - stick tmp file name to seq object

    my $seqobj = $self->_to_seqobj(%args);
    ($seqobj->{_fastaFH}, $seqobj->{_fastafile}) = tmpnam();
    # we need $fastafile below

    my $seqFH = Bio::SeqIO->newFh(-fh=>$seqobj->{_fastaFH}, -format=>"Fasta");
    print $seqFH $seqobj;


    # iterate through pwms
    my @PWMs;
    $self->reset();

    while (my $pwm = $self->next() ) {
	push @PWMs,$pwm;
    }

    # do the analysis

    my $hitlist = TFBS::SiteSet->new();

    foreach my $pwm (@PWMs)  {
	my $threshold = ($args{-threshold} or $pwm->{minscore});
	$hitlist->add_siteset($pwm->search_seq(-seqobj=>$seqobj,
					    -threshold =>$threshold ));
    }
    delete $seqobj->{_fastaFH};
    delete $seqobj->{_fastafile};
    return $hitlist;

    my ($self, $output) = @_;
    my ($optimal, $score_best, $in) = (undef, undef, 0);

    for (@$output) {
	if (/NumClust/) {
	    $in = 1;
	    next;
	}

	last if (/Tree Built/);

	if ($in) {
	    my (undef, $clusters, $score) = split(/\t/);

	    if ((not defined($score_best)) || $score < $score_best) {
		$score_best = $score;
		$optimal = $clusters;
	    }
	}
    }

    return $optimal;

    my ($self, %args) = @_;
    my $fh;

    for (keys(%args)) {
	die $error unless ($self->_check_opt($_, $args{$_}, $stamp_opt{$_}));
    }

    # Write matrices to temporary file
    if (not exists($args{-tf})) {
	$fh = new File::Temp( TEMPLATE => 'STAMP-XXXXX',
				 DIR => $args{-tempdir} || '/tmp',
				 SUFFIX => '.set');
	
	print $fh $_->STAMPprint() for (@{$self->{matrix_list}});
	$args{-tf} = $fh->filename();
    }
    # Set some default options
    $args{-tree} ||= "UPGMA";
    $args{-ma} ||= "IR";
    $args{-cc} ||= "PCC";
    $args{-align} ||= "SWU";


    # Make sure we find all files
    my $path;
    if ($args{-stampdir}) {
	$path = $args{-stampdir};
	die "Could not find STAMP at $path\n" if (not -e "$path/STAMP");
    } else {
	$path = (grep {-e "$_/STAMP"} split(/:+/, $ENV{PATH}))[0];
	$path || die "Could not find STAMP in path\n";
    }

    $args{-sd} ||= $path."/ScoreDists/JaspRand_".$args{-cc}."_".$args{-align}.".scores"; 
    die "No score distribution file found or not readable at '$args{-sd}'.\n Use -sd.\n" unless (-r $args{-sd});

    # Execute STAMP
    my $args = "";
    $args .= "$_ $args{$_} " for (keys(%args));
    my @output = `$path/STAMP -ch $args -out $fh`;
    
    # Get tree
    my $treeio = new Bio::TreeIO(-format => 'newick', -file =>  $fh->filename().".tree");
    my $tree = $treeio->next_tree;

    # Get FBP
    my $fbp = TFBS::Matrix::PFM->new(-matrixfile => $fh->filename()."FBP.txt");
    $fbp->{'filename'} = $fh->filename()."FBP.txt";

    print STDERR "::: $fbp->{'filename'}\n ";

    if (not $args{-noclean}) {
	my $deleted = unlink($fh->filename()."FBP.txt", $fh->filename().".tree");
	warn("Couldn't remove temporary files") if ($deleted != 2);
    }

    return ($fh,\@ output, $tree, $fbp);

    my ($self, %args) = @_;

    # DIRTY - stick tmp file name to seq object

    my $seqobj = $self->_to_seqobj(%args);
    ($seqobj->{_fastaFH}, $seqobj->{_fastafile}) = tmpnam();
    # we need $fastafile below

    my $outstream = Bio::SeqIO->new(-file=>">".$seqobj->{_fastafile}, -format=>"Fasta");
    my $subseqobj;
    if(my $subpart = $args{-subpart}) {
	my $subseq_start = $subpart->{-start};
	my $subseq_end = $subpart->{-end};
	unless($subseq_start and $subseq_end) {
	    $self->throw("Option -subpart missing suboption -relative_to, -start or -end");
	}
	$subseqobj = Bio::Seq->new(-seq => $seqobj->subseq($subseq_start, $subseq_end),
				    -id => $seqobj->id);
    }
    $outstream->write_seq($subseqobj or $seqobj);
    $outstream->close;

    # iterate through pwms
    my @PWMs;
    my $mxit = $self->Iterator();

    while (my $pwm = $mxit->next() ) {
	push @PWMs,$pwm;
    }

    # do the analysis

    my $hitlist = TFBS::SiteSet->new();

    foreach my $pwm (@PWMs)  {
	my $threshold = ($args{-threshold} or $pwm->{minscore});
	$hitlist->add_siteset($pwm->search_seq(-seqobj=>$seqobj,
					    -threshold =>$threshold,
					      -subpart=>$args{-subpart}));
    }
    delete $seqobj->{_fastaFH};
    unlink $seqobj->{_fastafile};
    delete $seqobj->{_fastafile};
    return $hitlist;

    # not OO - help function for search

    my $seqobj = shift;
    my $seqstring = uc($seqobj->seq());

    my @perlarray;
    foreach (qw(A C G T))  {
	my $seqtobits = $seqstring;
	eval "\$seqtobits =~ tr/$_/1/";  # curr. letter $_ to 1
	eval "\$seqtobits =~ tr/1/0/c";  # non-1s to 0
	push @perlarray, [split("", $seqtobits)];
    }
    return byte (\@perlarray);

    my ($self, %args) = @_;

    my $seq;
    if ($args{-file})  {    # not a Bio::Seq
	return Bio::SeqIO->new(-file => $args{-file},
			     -format => 'fasta',
			     -moltype => 'dna')->next_seq();
    }
    elsif ($args{-seqstring}
	   or $args{-seq})
    {   # I guess it's a string then
	return Bio::Seq->new(-seq  => ($args{-seqstring} or $args{-seq}),
			     -id => ($args{-seq_id} or "undefined"),
			     -moltype => 'dna');
    }
    elsif ($args{'-seqobj'} and ref($args{'-seqobj'}) =~ /Bio\:\:Seq/) {
	# do nothing (maybe check later)
	return $args{'-seqobj'};
    }
    #elsif (ref($format) =~ /Bio\:\:Seq/ and !defined $seq)  {
	# if only one parameter passed and it's a Bio::Seq
	#return $format;
    #}
    else  {
	$self->throw ("Wrong parametes passed to search method: ".%args);
    }


    # CONTINUE HERE TOMORROW

    my $self = shift;
    return $self->add_matrix(@_);

    my ($self, @matrices) = @_;
    foreach my $matrix (@matrices) {
	$self->throw("Argument to add_matrix_set not a TFBS::Matrix object")
	    unless $matrix->isa("TFBS::Matrix");
    }
    push @{$self->{matrix_list}}, @matrices;
    return $self;

    my ($self, @sets) = @_;
    foreach my $matrixset (@sets)  {
	$self->throw("Argument to add_matrix_set not a TFBS::Matrixset object")
	    unless $matrixset->isa("TFBS::MatrixSet");
	push @{$self->{matrix_list}}, @{$matrixset->{matrix_list}};
    }

    my ($self, %args) = @_;

    if ($self->size() <= 1) {
	warn("Can't cluster MatrixSet of size less than 2");
	return;
    }

    my ($fh, $output, $tree, $fbp) = $self->_run_STAMP(%args);

    # Find optimal cluster number
    my $optimal = $args{-optimal} || $self->_find_optimal($output);
    my $root = $tree->get_root_node();

    my @nodes = ($root);
    my @leaves;

    # Descend the tree until the optimal cluster number is reached
    while (scalar(@nodes) && (scalar(@nodes) + scalar(@leaves)) < $optimal) {
	my $node = pop @nodes;

	if ($node->is_Leaf()) {
	    push @leaves, $node;
	} else {	    
	    @nodes = sort {$a->height() <=> $b->height()} (@nodes, $node->each_Descendent());
	}
    } 

    # Build the clusters
    my @clusters;
    for (@leaves, @nodes) {
	my $cluster = $self->new();	
	$self->_build_cluster($cluster, $_);
	push @clusters, $cluster;
    }


    return ($tree, $optimal,\@ clusters);

    my ($self, %args) = @_;
    if ($self->size() == 0) {
	warn("Can't create FBP for MatrixSet of size 0");
	return;
    } elsif ($self->size() == 1) {
	return @{$self->{'matrix_list'}}[0];
    }

    my ($fh, $output, $tree, $fbp) = $self->_run_STAMP(%args);

    return $fbp;
}




1;

    my ($caller, @matrices) = shift;
    my $self = bless {matrix_list =>[]}, ref($caller) || $caller;
    $self->add_matrix(@matrices) if @matrices;
    return $self;

    my $class = shift;
    my %args = @_;
    my $self = bless {}, ref($class) || $class;

    if (defined $args{'-matrices'}) {
	$self->add_matrix( @{$args{'-matrices'}} ) if @{$args{'-matrices'}};
    } 
    
    if (defined $args{'-matrixfile'}) {
	my @matrices;

	open (FILE,  $args{-matrixfile})
	    or $self->throw("Could not open $args{-matrixfile}");

	while (<FILE>) {
	    /^\s*$/ && next;
	    if (/^>/) {
		
	    }
	}
	close(FILE);

	
    }

    return $self;

    my ($self) = @_;
    $self->warn("next: Deprecated method use Iterator instead.");
    if (my $next_matrix = shift (@{$self->{_iterator_list}})) {
	return $next_matrix;
    }
    else  {
	$self->reset;
	return undef;
    }

    my $self = shift;
    my (@lengths, @concat);
    my ($length, $i) = (-1, 0);

    # Concatenate to one big matrix
    for my $matrix (@{$self->{matrix_list}}) {
	$length += $matrix->length();
	push @lengths, $matrix->length();
	push @{$concat[$_]}, @{${$matrix->matrix()}[$_]} for (0..3);
    }

    # Schwartzian transform to get random permutation
    map { ( undef, $concat[0][$i], $concat[1][$i], $concat[2][$i],  $concat[3][$i] ) = @$_; $i++; } 
    sort { $a->[0] <=> $b->[0] } 
    map { [ rand(), $concat[0][$_], $concat[1][$_], $concat[2][$_], $concat[3][$_] ] } ( 0 .. $length );

    # Split it up again
    my $start = 0;
    for my $matrix (@{$self->{matrix_list}}) {
	my $length = shift(@lengths);
	my $end = $start + $length - 1;

	$matrix->matrix( [
			  [ @{$concat[0]}[$start..$end] ],
			  [ @{$concat[1]}[$start..$end] ],
			  [ @{$concat[2]}[$start..$end] ],
			  [ @{$concat[3]}[$start..$end] ]
			  ] 
			 );
	$start += $length;
    }

    my ($self, $id) = @_;

    my @list = grep { $_->ID() ne $id } @{$self->{matrix_list}};
    $self->{matrix_list} =\@ list;
}

my $error;

    my ($self) = @_;
    $self->warn("reset: Deprecated method use Iterator instead.");
    @{$self->{_iterator_list}} = @{$self->{matrix_list}};

    my ($self, %args) = @_;
    my $mxit = $self->Iterator();
    my $sitepairset = TFBS::SitePairSet->new;
    my $aln = TFBS::Matrix::_Alignment->new(%args);
    while (my $mx = $mxit->next) {
        my $singleset = $mx->search_aln(%args,
                                        -alignment_setup => $aln);
        $sitepairset->add_site_pair_set($singleset);
    }
    return $sitepairset;

    my ($self, %args) = @_;
    $self->_search(%args);

    scalar @{ $_[0]->{matrix_list} };

    my ($self) = @_;
    $self->warn("sort_by_name: Deprecated method use Iterator instead.");
    @{$self->{matrix_list}} = sort { uc($a->{name}) cmp uc ($b->{name}) }
                              @{$self->{matrix_list}};
    $self->reset();