~ubuntu-branches/ubuntu/lucid/pdl/lucid

« back to all changes in this revision

Viewing changes to Lib/Interpolate/Interpolate.pm

  • Committer: Bazaar Package Importer
  • Author(s): Ben Gertzfield
  • Date: 2002-04-08 18:47:16 UTC
  • Revision ID: james.westby@ubuntu.com-20020408184716-0hf64dc96kin3htp
Tags: upstream-2.3.2
ImportĀ upstreamĀ versionĀ 2.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
 
 
2
=head1 NAME
 
3
 
 
4
PDL::Interpolate - provide a consistent interface to the interpolation routines available in PDL
 
5
 
 
6
=head1 SYNOPSIS
 
7
 
 
8
 use PDL::Interpolate;
 
9
 
 
10
 my $i = new PDL::Interpolate( x => $x, y = $y );
 
11
 my $y = $i->interpolate( $xi ); 
 
12
 
 
13
=head1 DESCRIPTION
 
14
 
 
15
This module aims to provide a relatively-uniform interface
 
16
to the various interpolation methods available to PDL.
 
17
The idea is that a different interpolation scheme
 
18
can be used just by changing the C<new> call.
 
19
 
 
20
At present, PDL::Interpolate itself just provides
 
21
a somewhat-convoluted interface to the C<interpolate>
 
22
function of L<PDL::Primitive|PDL::Primitive/interpolate>. 
 
23
However, it is expected that derived classes,
 
24
such as 
 
25
L<PDL::Interpolate::Slatec|PDL::Interpolate::Slatec>,
 
26
will actually be used in real-world situations.
 
27
 
 
28
To use, create a PDL::Interpolate (or a derived class)
 
29
object, supplying it with its required attributes.
 
30
 
 
31
=head1 LIBRARIES
 
32
 
 
33
Currently, the avaliable classes are
 
34
 
 
35
=over 4
 
36
 
 
37
=item PDL::Interpolate
 
38
 
 
39
Provides an interface to the interpolation routines of PDL.
 
40
At present this is the linear interpolation routine
 
41
L<PDL::Primitive::interpol|PDL::Primitive/interpol>.
 
42
 
 
43
=item PDL::Interpolate::Slatec
 
44
 
 
45
The SLATEC library contains several approaches to interpolation:
 
46
piecewise cubic Hermite functions and B-splines.
 
47
At present, only the former method is available.
 
48
 
 
49
=back
 
50
 
 
51
It should be relatively easy to provide an interface to other
 
52
interpolation routines, such as those provided by the
 
53
Gnu Scientific Library (GSL).
 
54
 
 
55
=head1 ATTRIBUTES
 
56
 
 
57
The attributes (or options) of an object are as follows; 
 
58
derived classes may modify this list.
 
59
 
 
60
 Attribute  Flag  Description
 
61
 x          sgr   x positions of data
 
62
 y          sgr   function values at x positions
 
63
 bc         g     boundary conditions
 
64
 err        g     error flag
 
65
 type       g     type of interpolation
 
66
 
 
67
A flag of C<s> means that a user can set this attribute 
 
68
with the L<new|/new> or L<set|/set> methods,
 
69
a flag of C<g> means that the user can obtain the 
 
70
value of this attribute using L<get|/get>,
 
71
and a flag of C<r> means that the attribute is required
 
72
when an object is created (see the L<new|/new> method).
 
73
 
 
74
 Attribute  Default value
 
75
 bc         "none"
 
76
 type       "linear"
 
77
 
 
78
If a routine is sent an attribute it does not understand, then
 
79
it ignores that attribute, except for L<get|/get>, which
 
80
returns C<undef> for that value.
 
81
 
 
82
=head1 METHODS
 
83
 
 
84
The default methods are described below. However, defined classes
 
85
may extend them as they see fit, and add new methods.
 
86
 
 
87
Throughout this documentation, C<$x> and C<$y> refer to the function
 
88
to be interpolated whilst C<$xi> and C<$yi> are the interpolated values.
 
89
 
 
90
=head1 THREADING
 
91
 
 
92
The class will thread properly if the routines it calls do so.
 
93
See the SYNOPSIS section of L<PDL::Interpolate::Slatec>
 
94
(if available) for an example.
 
95
 
 
96
=cut
 
97
 
 
98
package PDL::Interpolate;
 
99
 
 
100
use Carp;
 
101
use strict;
 
102
 
 
103
####################################################################
 
104
 
 
105
## Public routines:
 
106
 
 
107
=head2 new
 
108
 
 
109
=for usage
 
110
 
 
111
 $obj = new PDL::Interpolate( x => $x, y => $y );
 
112
 
 
113
=for ref
 
114
 
 
115
Create a PDL::Interpolate object.
 
116
 
 
117
The required L<attributes|/attributes> are 
 
118
C<x> and C<y>.
 
119
At present the only available interpolation method 
 
120
is C<"linear"> - which just uses
 
121
L<PDL::Primitive::interpolate|PDL::Primitive::interpolate> - and
 
122
there are no options for boundary conditions, which is why
 
123
the C<type> and C<bc> attributes can not be changed.
 
124
 
 
125
=cut
 
126
 
 
127
# meaning of types:
 
128
#  required  - required, if this attr is changed, we need to re-initialise
 
129
#  settable  - can be changed with a new() or set() command
 
130
#  gettable  - can be read with a get() command
 
131
#
 
132
sub new {
 
133
    my $this  = shift;
 
134
    my $class = ref($this) || $this;
 
135
 
 
136
    # class structure
 
137
    my $self = { 
 
138
        attributes => {},
 
139
        values     => {},
 
140
        types      => { required => 0, settable => 0, gettable => 0 },
 
141
        flags      => { library => "PDL", status => 1, routine => "none", changed => 1 },
 
142
    };
 
143
 
 
144
    # make $self into an object
 
145
    bless $self, $class;
 
146
 
 
147
    # set up default attributes
 
148
    #
 
149
    $self->_add_attr( 
 
150
                      x    => { required => 1, settable => 1, gettable => 1 },
 
151
                      y    => { required => 1, settable => 1, gettable => 1 },
 
152
                      bc   => { gettable => 1 },
 
153
                      err  => { gettable => 1 },
 
154
                      type => { gettable => 1 }, 
 
155
                      );
 
156
    $self->_set_value( 
 
157
                       bc   => "none",
 
158
                       type => "linear",
 
159
                       );
 
160
 
 
161
    # set variables
 
162
    # - expect sub-classes to call this new with no variables, so $#_ == -1
 
163
    $self->set( @_ ) if ( @_ );
 
164
 
 
165
    # return the object
 
166
    return $self;
 
167
                                                                                
 
168
} # sub: new()
 
169
 
 
170
#####################################################################
 
171
 
 
172
# in _add_attr(), _change_attr() and _add_attr_type()
 
173
# we set flags->changed to 1 when something changes. It's
 
174
# a bit over the top to do this, as these should only be called when
 
175
# creating the object, when the changed flag should be set to 1 anyway
 
176
 
 
177
# add attributes to the object and sets value to undef
 
178
#
 
179
# supply a hash array, keys == variable name,
 
180
# values are a hash array with keys matching
 
181
# $self->{values}, which also gives the default value
 
182
# for the type
 
183
#
 
184
# this can only be used to create an attribute - 
 
185
# see _change_attr() to change an already exsiting attribute.
 
186
#
 
187
# the fields are set to the default values, then filled in with the supplied values
 
188
# any value that is unknown is ignored
 
189
#
 
190
sub _add_attr {
 
191
    my $self = shift;
 
192
    my %attrs = ( @_ );
 
193
 
 
194
    foreach my $attr ( keys %attrs ) {
 
195
        croak "ERROR: adding an attribute ($attr) which is already known.\n"
 
196
            if defined $self->{attributes}->{$attr};
 
197
 
 
198
        # set default values
 
199
        foreach my $type ( keys %{$self->{types}} ) {
 
200
            $self->{attributes}->{$attr}->{$type} = $self->{types}->{$type};
 
201
        }
 
202
 
 
203
        # change the values to those supplied
 
204
        foreach my $type ( keys %{$attrs{$attr}} ) {
 
205
            $self->{attributes}->{$attr}->{$type} = $attrs{$attr}->{$type}
 
206
                if exists $self->{types}->{$type};
 
207
        }
 
208
        # set value to undef
 
209
        $self->{values}->{$attr} = undef;
 
210
    }
 
211
    $self->{flags}->{changed} = 1;
 
212
 
 
213
} # sub: _add_attr()
 
214
 
 
215
# changes attributes of the object
 
216
 
217
# the given attributes MUST already exist
 
218
#
 
219
sub _change_attr {
 
220
    my $self = shift;
 
221
    my %attrs = ( @_ );
 
222
 
 
223
    foreach my $attr ( keys %attrs ) {
 
224
        croak "ERROR: changing an attribute ($attr) which is not known.\n"
 
225
            unless defined $self->{attributes}->{$attr};
 
226
 
 
227
        # change the values to those supplied
 
228
        foreach my $type ( keys %{$attrs{$attr}} ) {
 
229
            if ( exists $self->{types}->{$type} ) {
 
230
                $self->{attributes}->{$attr}->{$type} = $attrs{$attr}->{$type};
 
231
                $self->{flags}->{changed} = 1;
 
232
            }
 
233
        }
 
234
    }
 
235
} # sub: _change_attr()
 
236
 
 
237
# adds the given types to the allowed list, and
 
238
# updates all attributes to contain the default value
 
239
#
 
240
# Useful for sub-classes which add new types
 
241
#
 
242
sub _add_attr_type {
 
243
    my $self = shift;
 
244
    my %types = ( @_ );
 
245
 
 
246
    foreach my $type ( keys %types ) {
 
247
        croak "ERROR: adding type ($type) that is already known.\n"
 
248
            if exists $self->{types}->{$type};
 
249
        $self->{types}->{$type} = $types{$type};
 
250
 
 
251
        # loop through each attribute, adding this type
 
252
        foreach my $attr ( keys %{$self->{attributes}} ) {
 
253
            $self->{attributes}->{$attr}->{$type} = $types{$type};
 
254
        }
 
255
 
 
256
        $self->{flags}->{changed} = 1;
 
257
    }
 
258
} # sub: _add_attr_type()
 
259
 
 
260
####################################################################
 
261
 
 
262
# if an attribute has changed, check all required attributes
 
263
# still exist and re-initialise the object (for PDL::Interpolate
 
264
# this is a nop)
 
265
#
 
266
sub _check_attr {
 
267
    my $self = shift;
 
268
    return unless $self->{flags}->{changed};
 
269
 
 
270
    my @emsg;
 
271
    foreach my $name ( keys %{ $self->{attributes} } ) {
 
272
        if( $self->{attributes}->{$name}->{required} ) {
 
273
            push @emsg, $name unless defined($self->{values}->{$name});
 
274
        }
 
275
    }
 
276
    croak "ERROR - the following attributes must be supplied:\n [ @emsg ]\n"
 
277
        unless $#emsg == -1;
 
278
    
 
279
    $self->{flags}->{routine} = "none";
 
280
    $self->{flags}->{status} = 1;
 
281
    
 
282
    $self->_initialise;
 
283
    $self->{flags}->{new} = 0;
 
284
 
 
285
} # sub: check_attr()
 
286
 
 
287
####################################################################
 
288
#
 
289
# method to be over-ridden by sub-classes
 
290
 
 
291
# PDL::Interpolate needs no initialisation
 
292
#
 
293
sub _initialise {}
 
294
 
 
295
####################################################################
 
296
 
 
297
# a version of set that ignores the settable flag
 
298
# - for use by the class, not by the public
 
299
#
 
300
# it still ignores unknown attributes
 
301
#
 
302
sub _set_value {
 
303
    my $self = shift;
 
304
    my %attrs = ( @_ );
 
305
    
 
306
    foreach my $attr ( keys %attrs ) {
 
307
        if ( exists($self->{values}->{$attr}) ) {
 
308
            $self->{values}->{$attr} = $attrs{$attr};
 
309
            $self->{flags}->{changed} = 1;
 
310
        }
 
311
    }
 
312
 
 
313
} # sub: _set_value()
 
314
 
 
315
# a version of get that ignores the gettable flag
 
316
# - for use by the class, not by the public
 
317
#
 
318
# an unknown attribute returns an undef
 
319
#
 
320
sub _get_value {
 
321
    my $self = shift;
 
322
 
 
323
    my @ret;
 
324
    foreach my $name ( @_ ) {
 
325
        if ( exists $self->{values}->{$name} ) {
 
326
            push @ret, $self->{values}->{$name};
 
327
        } else {
 
328
            push @ret, undef;
 
329
        }
 
330
    }
 
331
 
 
332
    return wantarray ? @ret : $ret[0];
 
333
 
 
334
} # sub: _get_value()
 
335
 
 
336
####################################################################
 
337
 
 
338
=head2 set
 
339
 
 
340
=for usage
 
341
 
 
342
 my $nset = $obj->set( x = $newx, $y => $newy );
 
343
 
 
344
=for ref
 
345
 
 
346
Set attributes for a PDL::Interpolate object.
 
347
 
 
348
The return value gives the number of the supplied attributes
 
349
which were actually set. 
 
350
 
 
351
=cut
 
352
 
 
353
sub set {
 
354
    my $self = shift;
 
355
    my %vals = ( @_ );
 
356
 
 
357
    my $ctr = 0;
 
358
    foreach my $name ( keys %vals ) {
 
359
        if ( exists $self->{attributes}->{$name}->{settable} ) {
 
360
            $self->{values}->{$name} = $vals{$name};
 
361
            $ctr++;
 
362
        }
 
363
    }
 
364
 
 
365
    $self->{flags}->{changed} = 1 if $ctr;
 
366
    return $ctr;
 
367
 
 
368
} # sub: set()
 
369
 
 
370
####################################################################
 
371
 
 
372
=head2 get
 
373
 
 
374
=for usage
 
375
 
 
376
 my $x         = $obj->get( x );
 
377
 my ( $x, $y ) = $obj->get( qw( x y ) );
 
378
 
 
379
=for ref
 
380
 
 
381
Get attributes from a PDL::Interpolate object.
 
382
 
 
383
Given a list of attribute names, return a list of
 
384
their values; in scalar mode return a scalar value.
 
385
If the supplied list contains an unknown attribute,
 
386
C<get> returns a value of C<undef> for that
 
387
attribute.
 
388
 
 
389
=cut
 
390
 
 
391
sub get {
 
392
    my $self = shift;
 
393
 
 
394
    my @ret;
 
395
    foreach my $name ( @_ ) {
 
396
        if ( exists $self->{attributes}->{$name}->{gettable} ) {
 
397
            push @ret, $self->{values}->{$name};
 
398
        } else {
 
399
            push @ret, undef;
 
400
        }
 
401
    }
 
402
 
 
403
    return wantarray ? @ret : $ret[0];
 
404
 
 
405
} # sub: get()
 
406
 
 
407
####################################################################
 
408
 
 
409
=head2 interpolate
 
410
 
 
411
=for usage
 
412
 
 
413
 my $yi = $obj->interpolate( $xi );
 
414
 
 
415
=for ref
 
416
 
 
417
Returns the interpolated function at a given set of points.
 
418
 
 
419
A status value of -1, as returned by the C<status> method, 
 
420
means that some of the C<$xi> points lay outside the 
 
421
range of the data. The values for these points
 
422
were calculated using linear extrapolation.
 
423
 
 
424
=cut
 
425
 
 
426
sub interpolate {
 
427
    my $self = shift;
 
428
    my $xi   = shift;
 
429
 
 
430
    croak 'Usage: $obj->interpolate( $xi )' . "\n"
 
431
        unless defined $xi;
 
432
 
 
433
    # check everything is fine
 
434
    $self->_check_attr();
 
435
 
 
436
    # get values in one go
 
437
    my ( $x, $y ) = $self->_get_value( qw( x y ) );
 
438
 
 
439
    my ( $yi, $err ) = PDL::Primitive::interpolate( $xi, $x, $y );
 
440
 
 
441
    if ( any $err ) {
 
442
        $self->{flags}->{status} = -1;
 
443
    } else {
 
444
        $self->{flags}->{status} = 1;
 
445
    }
 
446
    $self->_set_value( err => $err );
 
447
 
 
448
    $self->{flags}->{routine} = "interpolate";
 
449
 
 
450
    return $yi;
 
451
}
 
452
 
 
453
####################################################################
 
454
#
 
455
# access to flags - have individual methods for these
 
456
 
 
457
=head2 status
 
458
 
 
459
=for usage
 
460
 
 
461
 my $status = $obj->status;
 
462
 
 
463
=for ref
 
464
 
 
465
Returns the status of a PDL::Interpolate object
 
466
 
 
467
Returns B<1> if everything is okay, B<0> if 
 
468
there has been a serious error since the last time
 
469
C<status> was called, and B<-1> if there
 
470
was a problem which was not serious.
 
471
In the latter case, C<$obj-E<gt>get("err")> may
 
472
provide more information, depending on the
 
473
particular class.
 
474
 
 
475
=cut
 
476
 
 
477
sub status { my $self = shift; return $self->{flags}->{status}; }
 
478
 
 
479
=head2 library
 
480
 
 
481
=for usage
 
482
 
 
483
 my $name = $obj->library;
 
484
 
 
485
=for ref
 
486
 
 
487
Returns the name of the library used by a PDL::Interpolate object
 
488
 
 
489
For PDL::Interpolate, the library name is C<"PDL">.
 
490
 
 
491
=cut
 
492
 
 
493
sub library { my $self = shift; return $self->{flags}->{library}; }
 
494
 
 
495
=head2 routine
 
496
 
 
497
=for usage
 
498
 
 
499
 my $name = $obj->routine;
 
500
 
 
501
=for ref
 
502
 
 
503
Returns the name of the last routine called by a PDL::Interpolate object.
 
504
 
 
505
For PDL::Interpolate, the only routine used is C<"interpolate">.
 
506
This will be more useful when calling derived classes,
 
507
in particular when trying to decode the values stored in the
 
508
C<err> attribute.
 
509
 
 
510
=cut
 
511
 
 
512
sub routine { my $self = shift; return $self->{flags}->{routine}; }
 
513
 
 
514
=head2 attributes
 
515
 
 
516
=for usage
 
517
 
 
518
 $obj->attributes;
 
519
 PDL::Interpolate::attributes;
 
520
 
 
521
=for ref
 
522
 
 
523
Print out the flags for the attributes of an object. 
 
524
Useful in case the documentation is just too opaque!
 
525
 
 
526
=for example
 
527
 
 
528
 PDL::Interpolate->attributes;
 
529
 Flags  Attribute
 
530
  SGR    x
 
531
  SGR    y
 
532
  G      err
 
533
  G      type
 
534
  G      bc                                                                      
 
535
 
 
536
=cut
 
537
 
 
538
# note, can be called with the class, rather than just
 
539
# an object
 
540
#
 
541
# to allow this, I've used a horrible hack - we actually
 
542
# create an object and then print out the attributes from that
 
543
# Ugh!
 
544
#
 
545
sub attributes { 
 
546
    my $self = shift;
 
547
 
 
548
    # ugh
 
549
    $self = $self->new unless ref($self);
 
550
 
 
551
    print "Flags  Attribute\n";
 
552
    while ( my ( $attr, $hashref ) = each %{$self->{attributes}} ) {
 
553
        my $flag = "";
 
554
        $flag .= "S" if $hashref->{settable};
 
555
        $flag .= "G" if $hashref->{gettable};
 
556
        $flag .= "R" if $hashref->{required};
 
557
        
 
558
        printf " %-3s    %s\n", $flag, $attr;
 
559
    }
 
560
    return;
 
561
}
 
562
 
 
563
####################################################################
 
564
 
 
565
=head1 AUTHOR
 
566
 
 
567
Copyright (C) 2000 Doug Burke (burke@ifa.hawaii.edu).
 
568
All rights reserved. There is no warranty. 
 
569
You are allowed to redistribute this software / documentation as 
 
570
described in the file COPYING in the PDL distribution.                                    
 
571
 
 
572
=head1 SEE ALSO
 
573
 
 
574
L<PDL>, perltoot(1).
 
575
 
 
576
=cut
 
577
 
 
578
####################################################################
 
579
# End with a true
 
580
1;