added c power to calc_fixedstep

git-svn-id: http://biopieces.googlecode.com/svn/trunk@323 74ccb610-7750-0410-82ae-013aeee3265d

diff --git a/code_perl/Maasha/C_bitarray.pm b/code_perl/Maasha/C_bitarray.pm
new file mode 100644 (file)
index 0000000..523cbae
--- /dev/null
+++ b/code_perl/Maasha/C_bitarray.pm
@@ -0,0 +1,193 @@
+package Maasha::C_bitarray;
+
+# Copyright (C) 2008 Martin A. Hansen.
+
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 2
+# of the License, or (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+# http://www.gnu.org/copyleft/gpl.html
+
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+# This modules make all functions for casting, filling, and scanning a C integer array.
+
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+
+use strict;
+
+use vars qw ( @ISA @EXPORT );
+
+@ISA = qw( Exporter );
+
+use Inline ( C => Config => DIRECTORY => $ENV{ 'BP_TMP' } );
+
+use Inline C => << 'END_C';
+
+#include <assert.h>
+
+
+void c_array_interval_fill( SV *array, int beg, int end, int score )
+{
+    /* Martin A. Hansen, November 2008. */
+
+    /* Incretments all values in a C array */
+    /* with a given score in a given interval. */
+
+    int i;
+    int len;
+    int size;
+
+    int *a = ( int * ) SvPV( array, len );
+
+    size = len / sizeof( int );
+
+    assert( beg >= 0 );
+    assert( end >= 0 );
+    assert( beg <= size );
+    assert( end < size );
+    assert( score > 0 );
+
+    for ( i = beg; i < end + 1; i++ ) {
+        a[ i ] += score;
+    }
+}
+
+
+void c_array_interval_scan( SV *array, int pos )
+{
+    /* Martin A. Hansen, November 2008. */
+
+    /* Locates the next interval of non-zero */
+    /* integers in a C array given a starting pos. */
+
+    int i;
+    int len;
+    int size;
+    int beg;
+    int end;
+
+    int *a = ( int * ) SvPV( array, len );
+
+    size = len / sizeof( int );
+
+    assert( pos >= 0 );
+    assert( pos <= size );
+
+    while ( pos < size && a[ pos ] == 0 ) { pos++; }
+
+    beg = pos;
+
+    while ( pos < size && a[ pos ] != 0 ) { pos++; }
+
+    end = pos - 1;
+
+    if ( beg > end )
+    {
+        beg = -1;
+        end = -1;
+    }
+
+    Inline_Stack_Vars;
+    Inline_Stack_Reset;
+
+    Inline_Stack_Push( sv_2mortal( newSViv( beg ) ) );
+    Inline_Stack_Push( sv_2mortal( newSViv( end ) ) );
+
+    Inline_Stack_Done;
+}
+
+
+void c_array_interval_get( SV *array, int beg, int end )
+{
+    /* Martin A. Hansen, November 2008. */
+
+    /* Get all values from an array in a given interval, */
+    /* and return these as a Perl list. */
+
+    int i;
+    int len;
+    int size;
+
+    int *a = ( int * ) SvPV( array, len );
+
+    size = len / sizeof( int );
+
+    assert( beg >= 0 );
+    assert( end >= 0 );
+    assert( beg <= size );
+    assert( end <= size );
+
+    Inline_Stack_Vars;
+    Inline_Stack_Reset;
+
+    for ( i = beg; i < end + 1; i++ ) {
+        Inline_Stack_Push( sv_2mortal( newSViv( a[ i ] ) ) );
+    }
+
+    Inline_Stack_Done;
+}
+
+
+void c_array_print( SV *array )
+{
+    /* Martin A. Hansen, November 2008. */
+
+    /* Outputs a C array to stdout. */
+
+    int i;
+    int len;
+
+    int *a = ( int * ) SvPV( array, len );
+
+    for ( i = 0; i < len / sizeof( int ); ++i ) {
+        printf( "%d\n", a[ i ] );
+    }
+}
+
+END_C
+ 
+
+sub c_array_init
+{
+    # Martin A. Hansen, November 2008.
+
+    # Initializes a zeroed C integer array using
+    # Perls vec function to create a bit
+    # vector.
+
+    my ( $size,   # number of elements in array
+         $bits,   # bit size
+       ) = @_;
+
+    # Returns a bit vector
+
+    my ( $vec );
+
+    $vec  = '';
+
+    vec( $vec, $size - 1, $bits ) = 0;
+
+    return $vec;
+}
+
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+1;

diff --git a/code_perl/Maasha/UCSC/Wiggle.pm b/code_perl/Maasha/UCSC/Wiggle.pm
index 1fbd926c3f5f74238b5a359a1f2adcec93a49172..40bdc40dec2142c872dc90228d421d7c05fc1927 100644 (file)
--- a/code_perl/Maasha/UCSC/Wiggle.pm
+++ b/code_perl/Maasha/UCSC/Wiggle.pm
@@ -35,6 +35,7 @@ use strict;
 use Data::Dumper;
 use Maasha::Common;
 use Maasha::Filesys;
+use Maasha::C_bitarray;
 use Maasha::Calc;
 
 use vars qw( @ISA @EXPORT_OK );
@@ -48,7 +49,9 @@ require Exporter;
 
 
 use constant {
-    chrom       => 0,   # BED field names
+    BITS        => 32,            # Number of bits in an integer
+    SEQ_MAX     => 250_000_000,   # Maximum sequence size
+    chrom       => 0,             # BED field names
     chromStart  => 1,
     chromEnd    => 2,
     name        => 3,
@@ -213,9 +216,9 @@ sub fixedstep_calc
 
     $beg = 0;
 
-    while ( ( $beg, $end ) = fixedstep_scan( $array, $beg ) and $beg )
+    while ( ( $beg, $end ) = Maasha::C_bitarray::c_array_interval_scan( $array, $beg ) and $beg != -1 )
     {
-        @block = @{ $array }[ $beg .. $end - 1 ];
+        @block = Maasha::C_bitarray::c_array_interval_get( $array, $beg, $end );
 
         map { $_ = sprintf "%.4f", Maasha::Calc::log10( $_ ) } @block if $use_log10;
 
@@ -223,11 +226,13 @@ sub fixedstep_calc
 
         fixedstep_entry_put( \@block, $fh_out );
 
-        $beg = $end;
+        $beg = $end + 1;
     }
 
     close $fh_out;
 
+    undef $array;
+
     return $fixedstep_file;
 }
 
@@ -246,9 +251,11 @@ sub fixedstep_calc_array
          $use_score,   # flag indicating that the score field should be used - OPTIONAL
        ) = @_;
 
-    # Returns arrayref.
+    # Returns a bitarray.
+
+    my ( $bed, $clones, $vec );
 
-    my ( $bed, $clones, @array );
+    $vec = Maasha::C_bitarray::c_array_init( SEQ_MAX, BITS );
 
     while ( $bed = Maasha::UCSC::BED::bed_entry_get( $fh_in, 5 ) )
     {
@@ -260,10 +267,10 @@ sub fixedstep_calc_array
             $clones = 1;
         }
 
-        map { $array[ $_ ] += $clones } $bed->[ chromStart ] .. $bed->[ chromEnd ] - 1;
+        Maasha::C_bitarray::c_array_interval_fill( $vec, $bed->[ chromStart ], $bed->[ chromEnd ] - 1, $clones );
     }
 
-    return wantarray ? @array : \@array;
+    return $vec;
 }