From 2f764d2bab315006d4f2d96060152ddb711b0c5a Mon Sep 17 00:00:00 2001 From: martinahansen Date: Thu, 27 Nov 2008 02:57:52 +0000 Subject: [PATCH] added c power to calc_fixedstep git-svn-id: http://biopieces.googlecode.com/svn/trunk@323 74ccb610-7750-0410-82ae-013aeee3265d --- code_perl/Maasha/C_bitarray.pm | 193 ++++++++++++++++++++++++++++++++ code_perl/Maasha/UCSC/Wiggle.pm | 23 ++-- 2 files changed, 208 insertions(+), 8 deletions(-) create mode 100644 code_perl/Maasha/C_bitarray.pm diff --git a/code_perl/Maasha/C_bitarray.pm b/code_perl/Maasha/C_bitarray.pm new file mode 100644 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 + + +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 1fbd926..40bdc40 100644 --- 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; } -- 2.39.5