pwiz::math Namespace Reference

Namespaces
namespace	MatchedFilter
namespace	types

Classes
class	Random
class	LinearLeastSquares< LinearLeastSquaresType_LU >
class	LinearLeastSquares< LinearLeastSquaresType_QR >
class	LinearSolver< LinearSolverType_LU >
class	LinearSolver< LinearSolverType_QR >
struct	OrderedPair
class	OrderedPairContainerRef
	wrapper class for accessing contiguous data as a container of OrderedPairs; note that it does not own the underlying data More...
class	Parabola
class	Stats

Functions
PWIZ_API_DECL double	erf (double x)
	real error function; calls gcc-provided erf, complex version (below) on msvc
PWIZ_API_DECL std::complex < double >	erf (const std::complex< double > &z)
	complex error function
std::complex< double >	erf_series2 (const std::complex< double > &z)
	series implementation for testing
template<class T >
void	TransposeMultiply (const ublas::vector< T > &vector, ublas::matrix< T > &result, size_t size)
template<class T >
void	HouseholderCornerSubstraction (ublas::matrix< T > &LeftLarge, const ublas::matrix< T > &RightSmall)
template<class T >
void	HouseholderQR (const ublas::matrix< T > &M, ublas::matrix< T > &Q, ublas::matrix< T > &R)
std::ostream &	operator<< (std::ostream &os, const OrderedPair &p)
std::istream &	operator>> (std::istream &is, OrderedPair &p)
bool	operator== (const OrderedPair &a, const OrderedPair &b)
bool	operator!= (const OrderedPair &a, const OrderedPair &b)
PWIZ_API_DECL std::ostream &	operator<< (std::ostream &os, const Parabola &p)
template<class matrix_type , class vector_type >
void	Reflector (const vector_type &x, matrix_type &F)
template<class matrix_type , class vector_type >
matrix_type	Reflector (const vector_type &x)
template<class matrix_type >
void	qr (const matrix_type &A, matrix_type &Q, matrix_type &R)

---

Function Documentation

PWIZ_API_DECL double pwiz::math::erf

(

double

x

)

real error function; calls gcc-provided erf, complex version (below) on msvc

Referenced by test_real(), test_real_wrapper(), and test_series().

PWIZ_API_DECL std::complex<double> pwiz::math::erf

(

const std::complex< double > &

z

)

complex error function

std::complex<double> pwiz::math::erf_series2

(

const std::complex< double > &

z

)

series implementation for testing

Referenced by test_series().

template<class T >

void pwiz::math::TransposeMultiply	(	const ublas::vector< T > &	vector,
		ublas::matrix< T > &	result,
		size_t	size
	)

Definition at line 39 of file HouseholderQR.hpp.

Referenced by HouseholderQR().

{
  result.resize (size,size);
  result.clear ();
  for(unsigned int row=0; row< vector.size(); ++row)
    {
      for(unsigned int col=0; col < vector.size(); ++col)
        result(row,col) = vector(col) * vector(row);

    }
}

template<class T >

void pwiz::math::HouseholderCornerSubstraction	(	ublas::matrix< T > &	LeftLarge,
		const ublas::matrix< T > &	RightSmall
	)

Definition at line 54 of file HouseholderQR.hpp.

Referenced by HouseholderQR().

{
  using namespace boost::numeric::ublas;
  using namespace std; 
  if( 
     !( 
       (LeftLarge.size1() >= RightSmall.size1())
       && (LeftLarge.size2() >= RightSmall.size2())
        ) 
      )
    {
      cerr << "invalid matrix dimensions" << endl;
      return;
    }  

  size_t row_offset = LeftLarge.size2() - RightSmall.size2();
  size_t col_offset = LeftLarge.size1() - RightSmall.size1();

  for(unsigned int row = 0; row < RightSmall.size2(); ++row )
    for(unsigned int col = 0; col < RightSmall.size1(); ++col )
      LeftLarge(col_offset+col,row_offset+row) -= RightSmall(col,row);
}

template<class T >

void pwiz::math::HouseholderQR	(	const ublas::matrix< T > &	M,
		ublas::matrix< T > &	Q,
		ublas::matrix< T > &	R
	)

Definition at line 79 of file HouseholderQR.hpp.

References H, HouseholderCornerSubstraction(), and TransposeMultiply().

Referenced by main().

{
  using namespace boost::numeric::ublas;
  using namespace std;  

  if( 
     !( 
       (M.size1() == M.size2())
        ) 
      )
    {
      cerr << "invalid matrix dimensions" << endl;
      return;
    }
  size_t size = M.size1();

  // init Matrices
  matrix<T> H, HTemp;
  HTemp = identity_matrix<T>(size);
  Q = identity_matrix<T>(size);
  R = M;

  // find Householder reflection matrices
  for(unsigned int col = 0; col < size-1; ++col)
    {
      // create X vector
      ublas::vector<T> RRowView = column(R,col);      
      vector_range< ublas::vector<T> > X2 (RRowView, range (col, size));
      ublas::vector<T> X = X2;

      // X -> U~
      if(X(0) >= 0)
        X(0) += norm_2(X);
      else
        X(0) += -1*norm_2(X);      

      HTemp.resize(X.size(),X.size(),true);

      TransposeMultiply(X, HTemp, X.size());

      // HTemp = the 2UUt part of H 
      HTemp *= ( 2 / inner_prod(X,X) );

      // H = I - 2UUt
      H = identity_matrix<T>(size);
      HouseholderCornerSubstraction(H,HTemp);

      // add H to Q and R
      Q = prod(Q,H);
      R = prod(H,R);
    }
}

std::ostream& pwiz::math::operator<< ( std::ostream &  os, const OrderedPair &  p  )

inline

Definition at line 49 of file OrderedPair.hpp.

References pwiz::math::OrderedPair::x, and pwiz::math::OrderedPair::y.

{
    os << "(" << p.x << "," << p.y << ")";
    return os;
}

std::istream& pwiz::math::operator>> ( std::istream &  is, OrderedPair &  p  )

inline

Definition at line 56 of file OrderedPair.hpp.

References pwiz::math::OrderedPair::x, and pwiz::math::OrderedPair::y.

{
    char open='\0', comma='\0', close='\0';
    is >>  open >> p.x >> comma >> p.y >> close;
    if (!is) return is;
    if (open!='(' || comma!=',' || close!=')')
        throw std::runtime_error("[OrderedPair::operator>>] Unexpected input.");
    return is;
}

bool pwiz::math::operator== ( const OrderedPair &  a, const OrderedPair &  b  )

inline

Definition at line 67 of file OrderedPair.hpp.

References pwiz::math::OrderedPair::x, and pwiz::math::OrderedPair::y.

{
    return a.x==b.x && a.y==b.y;
}

bool pwiz::math::operator!= ( const OrderedPair &  a, const OrderedPair &  b  )

inline

Definition at line 73 of file OrderedPair.hpp.

{
    return !(a == b);
}

PWIZ_API_DECL std::ostream& pwiz::math::operator<<	(	std::ostream &	os,
		const Parabola &	p
	)

template<class matrix_type , class vector_type >

void pwiz::math::Reflector	(	const vector_type &	x,
		matrix_type &	F
	)

Definition at line 39 of file qr.hpp.

References x.

Referenced by qr(), and testReflector().

{
    using namespace boost::numeric::ublas;

    typedef typename matrix_type::value_type value_type;

    unit_vector<value_type> e1(x.size(), 0);

    //v_k = -sgn( x(1) ) * inner_prod(x) * e1 + x;
    double x_2 = norm_2(x);
    boost::numeric::ublas::vector<value_type>
        v_k((x(0) >= 0 ? x_2 : -1 * x_2) * e1 + x);

    //v_k = v_k / norm_2(v_k);
    double norm_vk = norm_2(v_k);
    if (norm_vk != 0)
        v_k /= norm_2(v_k);
    
    // F = A(k:m,k:n) - 2 * outer_prod(v_k, v_k) * A(k:m,k:n)
    identity_matrix<value_type> eye(v_k.size());
    F = matrix_type(v_k.size(), v_k.size());
    
    F = eye - 2. * outer_prod(v_k, v_k);
}

template<class matrix_type , class vector_type >

matrix_type pwiz::math::Reflector

(

const vector_type &

x

)

Definition at line 69 of file qr.hpp.

References F, and x.

{
    using namespace boost::numeric::ublas;

    matrix_type F(x.size(), x.size());

    Reflector<matrix_type, vector_type>(x, F);

    return F;
}

template<class matrix_type >

void pwiz::math::qr	(	const matrix_type &	A,
		matrix_type &	Q,
		matrix_type &	R
	)

Definition at line 81 of file qr.hpp.

References A, F, ident(), k, Reflector(), and x.

Referenced by pwiz::math::LinearSolver< LinearSolverType_QR >::solve().

{
    using namespace boost::numeric::ublas;

    typedef typename matrix_type::size_type size_type;
    typedef typename matrix_type::value_type value_type;

    // TODO resize Q and R to match the needed size.
    int m=A.size1();
    int n=A.size2();

    identity_matrix<value_type> ident(m);
    if (Q.size1() != ident.size1() || Q.size2() != ident.size2())
        Q = matrix_type(m, m);
    Q.assign(ident);

    R.clear();
    R = A;

    for (size_type k=0; k< R.size1() && k<R.size2(); k++)
    {
        slice s1(k, 1, m - k);
        slice s2(k, 0, m - k);
        unit_vector<value_type> e1(m - k, 0);

        // x = A(k:m, k);
        matrix_vector_slice<matrix_type> x(R, s1, s2);
        matrix_type F(x.size(), x.size());
        
        Reflector(x, F);

        matrix_type temp = subrange(R, k, m, k, n);
        //F = prod(F, temp);
        subrange(R, k, m, k, n) = prod(F, temp);

        // <<---------------------------------------------->>
        // forming Q
        identity_matrix<value_type> iqk(A.size1());
        matrix_type Q_k(iqk);
        
        subrange(Q_k, Q_k.size1() - F.size1(), Q_k.size1(),
                 Q_k.size2() - F.size2(), Q_k.size2()) = F;

        Q = prod(Q, Q_k);
    }
}