CHAPTER 2

Math Utilities

 

Several math utilities, constants, pointers, data structures are defined in MAC and other sub-modules of CPL.  They do not depend on any other parts of the libraries.

 

2.1 utility.h and cpl.mac.utility

 

In C++, a utility.h is defined and subsequently shadowed into Python.

 

2.1.1 In C++

 

The following functions are defined under the namespace of mac.

 

double rand_fractional(const unsigned int seed = 0);

rand_fractional returns a random fractional number between 0 and 1.  It may accept an unsigned int as a seed.

 

template <class T> bool even(const T i);

template <class T> bool  odd(const T i);

The template functions even and odd test whether an integer value is even or odd, respectively.

 

template <class T> bool feven(const T,

                              const T = TINY);

template <class T> bool  fodd(const T,

                              const T = TINY);

The template functions feven and fodd test whether a float point value is close to even or odd, respectively.  They may accept a small number as tolerance.  The default tolerance is TINY defined to 1e-6.

 

template <class T> int nearest(const T x);

The template function nearest returns an int that is closest to the given value.

 

template <class T> std::complex<int> nearestComplex(

   const std::complex<T>);

Returns a complex<int> that is closest to the given complex.

 

template <class T> bool approx(const T& x,

                               const T& y = 0,

                               const double t = TINY);

Returns true, if the first two objects are approximately equal to each other.  The optional double is a tolerance.  ¡®Approximately equal¡¯ is defined as (x ¨C y) * (x ¨C y) < t * t, where operator ¨C must be defined for type T to return the difference.  The operator * must be defined for type T to perform multiplication or to return dot product.  The return type shall be scalar such as double.  This definition does not work with complex<T>, since the operator * does not meet the requirement.  Use the function approxComplex().

 

template <class T> bool approxComplex(

   const T& x,

   const T& y = 0,

   const double t = TINY);

Returns true, if the first two complex objects are approximately equal to each other.  The optional double is a tolerance.  Operator ¨C and function norm() must be defined for type T.

 

template <class X, class Y>

Y ramp(const X&, const X&, const Y&,

                 const X&, const Y&);

The template function ramp(x, x1, y1, x2, y2) computes a linear interpolation at x on the straight line from (x1, y1) to (x2, y2) or on the extension of the straight line.  Obviously, x1 and x2 must be different.  If they are the same, an assertion error will occur.  All x¡¯s may have a different type from the y¡¯s.

 

template <class T> T pythag(const T&, const T&);

The template function pythag(a, b) computes  without destructive underflow and overflow.  a and b obviously cannot be of integer types.  If they are, first the compilation may not pass; second, even if it passes OK, the result may not be what the user intends.

 

The following macros are defined, but they are not under any namespace.

 

degreeToRadian

(M_PI / 180)

radianToDegree

(180 / M_PI)

TINY

1e-6

TEENY

1e-15

TEENYTINY

1e-300

SIGN(a,b)

((b) >= 0.0 ? fabs(a) : -fabs(a))

 

Table 2.1.1.0.1 Macros defined in utility.h.

 

template <class T = double> class lessComplex

A binary functor that makes complex less-comparable, so that complex can be used as a key in set or map.  If the real component is equivalent, the imaginary component defines ¡®less than¡¯, otherwise, the real component defines ¡®less than¡¯.

 

2.1.2 In Python

 

All functions and macros above except SIGN(a,b) are also available in Python as members of the module cpl.mac.utility.

 

rand_fractional(seed = 0)

Returns a random fractional number.  It may take an unsigned integer as a seed for random number generation.

 

even(i)

odd(i)

Test whether the given integer i is even or odd, respectively.

 

feven(f, t = cpl.mac.utility.TINY)

fodd(f, t = cpl.mac.utility.TINY)

Test whether the given floating number f is an even or odd integer, respectively.  The testing is done within the optional tolerance t.

 

nearest(f)

Returns the nearest integer of the given floating number f.

 

nearestComplex(c)

Returns the nearest complex of integer to the given complex c.

 

approx(x, y = 0, t = cpl.mac.utility.TINY)

Returns true, if x and y are approximately equal to each other within the tolerance of t.  x and y cannot be complex.

 

approxComplex(x, y = 0, t = cpl.mac.utility.TINY)

Returns true, if two complexes x and y are approximately equal to each other within the tolerance of t.  x and y cannot be scalars.

 

ramp(x, x1, y1, x2, y2)

pythag(a, b)

See functions ramp() and pythag() in 2.1.1 for detail.  All variables and the return value are double.

 

The constants listed in Table 2.1.1.0.1 can be accessed by prefix cpl.mac.utility to the constant name.

 

2.2 Accessing std data structures from CPL

 

The standard template library (STL) of C++ contains a set of very useful data structures such as std::vector<T>.  They are commonly used in C++ programs.  CCL uses them extensively.  In order to wrap the CCL classes into Python, communication between the basic data structures of STL and those in Python is essential.  Thanks to some available typemaps of SWIG (http://www.swig.org), such task becomes trivial.

 

2.2.1 std::vector<T>

 

std::vector<T> is instantiated with T being the built-in types double, float, long int, int, short int, long unsigned int, unsigned int, short unsigned int, and bool.  They can be accessed by using the following class names:

 

cpl.std.vectorLDC

cpl.std.vectorDCX

cpl.std.vectorFCX

cpl.std.vectorLDB

cpl.std.vectorDBL

cpl.std.vectorFLT

cpl.std.vectorLIN

cpl.std.vectorINT

cpl.std.vectorSIN

cpl.std.vectorLUI

cpl.std.vectorUIN

cpl.std.vectorSUI

cpl.std.vectorBOL

cpl.std.vectorCHA

cpl.std.vectorSTR

 

All of them have the same API.  In this section, cpl.std.vectorTYP represents any of them above.

 

vectorTYP(n = 0)

Constructs a vector vectorTYP of size n with all 0 elements.

 

The rest of the member functions and operators resemble those of a Python list.

 

When T is the type defined in CCL or MAC, the vector of such type is located in the same module as the class T would locate, for example, cpl.ccl.symm.MillerIndices and cpl.ccl.symm.vectorH.

 

2.2.2 std::map<K, T>

to be developed

 

2.3 Accessing SWIG pointer library from CPL

 

The SWIG pointer library offers a means to communicate between a C array to and from Python sequence such as list, and ways to access directly to specific memory in C from Python.  The SWIG pointer library is wrapped in CPL as module cpl.ptr.

 

2.3.1 Pointer class

 

A pointer class acts like a pointer, but it is a class.  See the SWIG manual for documentation at http://www.swig.org/Doc1.3/Library.html.  They are named pointerTYP in module cpl.ptr, where TYP can be DBL, FLT, LIN, INT, SIN, LUI, UIN, SUI, BOL, and CHA (see Table 1.3.2.0.1).

 

pointerTYP()

Constructs a pointer object that points to a new object of type TYP.  The newly-constructed pointer (as p in the rest of this section) can be passed to C functions that expect such pointer.

 

p.assign(a)

Assigns the given value a of type TYP to the object pointed by p.

 

p.value()

Returns the value that is pointed by p.

 

p.cast()

Returns a bare pointer (not the pointer object) to the value pointed by p.

 

pointerTYP.frompointer(p)

Creates a new pointer object that points to the existing object that is pointed by p.  Here p can be either a bare pointer or a pointer object, but it must point to type TYP.

 

#!/usr/bin/python

 

import ptr

 

a = 7

p = ptr.pointerINT()              # p is a new pointer object.

p.assign(a)                       # The object pointed by p is assigned by a.

assert p.value() == a

 

t = p.cast()                      # t is bare pointer to value a.

q = ptr.pointerINT.frompointer(t) # q is a new pointer object different than p,

                                  # but q points to the same object as pointed

                                  # by t.

assert p != q

assert p.value() == q.value()

 

q.assign(-a)                      # The object pointed by q is assigned by a

                                  # new value.  So is the object pointed by p.

assert p != q

assert p.value() == q.value()

 

Listing 2.3.1.0.1 Testing pointer class.

 

2.3.2 Array class

 

An object of array class wraps a C pointer that points to an array, and be passed to C functions that expect an array.  They are named arrayTYP in module cpl.ptr, where TYP can be DBL, FLT, LIN, INT, SIN, LUI, UIN, SUI, BOL, and CHA (see Table 1.3.2.0.1).

 

arrayTYP(n)

Constructs an array object of size n, called a in the rest of this section.  The initial values in the array are arbitrary.

 

a[i]

Can be assigned by a value of type TYP, and can return the value of the element with the index i.

 

a.cast()

Returns a bare pointer to the array.

 

arrayTYP.frompointer(p)

Creates a new array object that wraps the existing array pointed by p.  Here p can be either a bare pointer, a pointer object or an array object, but it must point to type TYP.

 

#!/usr/bin/python

 

import ptr

 

n = 7

a = ptr.arrayINT(n)             # a is a new array object with n elements.

for i in range(n):

   a[i] = i * i                 # Assign values to the array.

for i in range(n):

   assert a[i] == i * i         # Get values from the array.

 

p = a.cast()                    # p is a bare pointer to the array.

b = ptr.arrayINT.frompointer(p) # b is a new name for the array a.

assert a != b

for i in range(n):

   assert a[i] == b[i]

 

for i in range(n):

   b[i] += i                    # The array b is changed.  So is array a.

assert a != b

for i in range(n):

   assert a[i] == b[i]

 

Listing 2.3.1.0.1 Testing array class.

 

2.3.3 Pointer functions

 

SWIG pointer library also contains various functions related to pointers to type TYP, where TYP can be DBL, FLT, LIN, INT, SIN, LUI, UIN, SUI, BOL, and CHA (see Table 1.3.2.0.1).  These types and pointer to these types are available.  They are named pTYP and ppTYP, respectively.

 

new_pTYP()

Creates a new object of type TYP and returns a pointer to the new object.

 

copy_pTYP(a)

Creates a new object of type TYP and returns a pointer to the new object.  The value of the new object is copied from a of type TYP.

 

delete_pTYP(p)

Deletes the object of type TYP pointed by p.

 

pTYP_assign(p, a)

Assigns value a of type TYP to the object pointed by p.

 

pTYP_value(p)

Returns the value of the object pointed by p.

 

#!/usr/bin/python

 

import ptr

 

a = 7

p = ptr.new_pINT()            # p is a new pointer object.

ptr.pINT_assign(p, a)         # The object pointed by p is assigned by a.

assert ptr.pINT_value(p) == a # Get value of object pointed by p.

 

q = ptr.copy_pINT(a)          # q is a new pointer object.  The value pointed

                              # by q is copied from a.

assert p != q

assert ptr.pINT_value(p) == ptr.pINT_value(q)

 

ptr.delete_pINT(p)            # The object pointed by p is deleted, but that

                              # pointed by q is still intact.

assert ptr.pINT_value(q) == a

 

ptr.delete_pINT(q)            # The object pointed by q is also deleted.

 

Listing 2.3.3.0.1 Testing pointer functions.

 

See also Listing 6.1.2.0.1 for an example for ppDBL.

 

2.3.4 Array functions

 

More functions are available to work with pointers that point to arrays.  Up to 3-dimensional arrays are available (aTYP, aaTYP, and aaaTYP), but these functions are not very convenient for multi-dimensional array operation; they are also error prone.  Available types are DBL, FLT, LIN, INT, SIN, LUI, UIN, SUI, BOL, and CHA (see Table 1.3.2.0.1).

 

new_aTYP(n)

Creates a new n-element array of type TYP.

 

delete_aTYP(p)

Deletes the array of type TYP pointed by p.

 

aTYP_setitem(p, i, a)

Assigns the value a of type TYP to the element of the array pointed by p at the 0-offset index i.

 

aTYP_getitem(p, i)

Returns the value of the element in the array pointed by p at the 0-offset index i.

 

#!/usr/bin/python

 

import ptr

 

m = 10

n = 7

a = ptr.new_aaINT(m)

         # a is a pointer to a new n-row 2-dimensional array.

for i in range(m):

   ptr.aaINT_setitem(a, i, ptr.new_aINT(n))

         # Each row of array a is a pointer to a 1-dimensional array.

   for j in range(n):

      ptr.aINT_setitem(ptr.aaINT_getitem(a, i), j, i + j)

         # Values are assigned to each row.

 

for i in range(m):

   for j in range(n):

      assert ptr.aINT_getitem(ptr.aaINT_getitem(a, i), j) == i + j

 

Listing 2.3.4.0.1 Testing array functions that construct a 2-dimensional array.

 

2.3.5 Conversion between C array and Python sequence

 

Module cpl.mathematics uses SWIG pointer library to define several functions to convert back and forth between a C array and a Python list, which make lots of common tasks much easier than using the pointer library directly.

 

arrayCreate(n = 0, type = "double", init = 0)

Creates and returns a new n-element array object of type type filled by the initial value init.  All types listed in Table 1.3.2.0.1 are valid.  type is a string from both sides of the table.  The array object could then be passed to C functions that require a pointer to array.

 

sequenceToArray(s, type = ¡¯double¡¯)

Returns an array object created by cpl.ptr.arrayTYP that contains the items of the given sequence s.  The sequence can be either a list or a tuple as long as the sequence supports operator [] for element access.  The data type may be specified by type.

 

arrayToList(a, n)

Returns a list that contains the items of the n-element array object a.  n can be less than the real size of array a.  If so, the first n elements are copied to the list.  This function works for any type of array element.

 

pointerToList(p, n, type = ¡¯double¡¯)

Returns a list that contains the items of the n-element array pointed by p.  n can be less than the real size of array a.  If so, the first n elements are copied to the list.  The type of the array element needs to be specified.

 

#!/usr/bin/python

 

import cpl

 

n = 7

a = cpl.mathematics.arrayCreate(n, 'int')     # a is a new array.

for i in range(n): a[i] = i * i               # Assign values to array a.

l = cpl.mathematics.arrayToList(a, n)         # l is a new list that contains

                                              # the same data as array a.

for i in range(n): assert a[i] == l[i]

 

b = cpl.mathematics.sequenceToArray(l, 'int') # b is new array that contains

                                              # the same data as array a.

assert a != b

for i in range(n): assert a[i] == b[i]

 

for i in range(n): b[i] += n                  # b is changed, but a is not.

assert a != b

for i in range(n): assert a[i] - b[i] == -n

 

k = cpl.mathematics.pointerToList(a.cast(), n, 'int')

   # k is a new list that contains the same data as array a.  k is constructed

   # from a pointer to array a.

assert k == l

 

Listing 2.3.5.0.1 Testing conversion between C array and Python sequence.