matvecop.hpp

模糊聚類分析源碼。包含教學文件

template <class D,class M1> inline 
void matrix_scale(M1& A, const D s)
  {
  matrix_transform
    (
    A,
    set_to<typename M1::value_type,typename M1::value_type>(),
    bind1st(multiplies<typename M1::value_type>(),s),
    A
    ); 
  }

template <class M1,class D,class M2> inline 
void matrix_set_scaled(M1& A, const D s, const M2& B)
  {
  matrix_transform
    (
    A,
    set_to<typename M1::value_type,typename M2::value_type>(),
    bind1st(multiplies<typename M2::value_type>(),s),
    B
    ); 
  }

template <class M1,class D,class M2> inline 
void matrix_inc_scaled(M1& A,const D s, const M2& B)
  { 
  matrix_transform
    (
    A,
    incr_by<typename M1::value_type,typename M2::value_type>(),
    bind1st(multiplies<typename M2::value_type>(),s),
    B
    ); 
  }

template <class M1,class D,class M2> inline 
void matrix_dec_scaled(M1& A,const D s, const M2& B)
  { 
  matrix_transform
    (
    A,
    decr_by<typename M1::value_type,D>(),
    bind1st(multiplies<typename M2::value_type>(),s),
    B
    ); 
  }

template <class M1,class M2> inline 
void matrix_div(M1& A,const M2& B)
  {
  matrix_transform
    (
    A,
    div_by<typename M1::value_type,typename M2::value_type>(),
    identity<typename M2::value_type>(),
    B
    ); 
  }

template <class M1,class D1,class M2,class D2> inline 
void matrix_set_combi(M1& A,const D1 s,const M2& B,const D2 t,const M2& C)
  { 
  matrix_transform
    (
    A,
    set_to<typename M1::value_type,typename M2::value_type>(),
    B,
    compose
      (
      plus<typename M2::value_type>(),
      bind2nd(multiplies<typename M2::value_type>(),s),
      bind2nd(multiplies<typename M2::value_type>(),t) 
      ),
    C
    ); 
  }

template <class M1,class D1,class M2,class D2> inline 
void matrix_inc_combi(M1& A,const D1 s,const M2& B,const D2 t,const M2& C)
  { 
  matrix_transform
    (
    A,
    incr_by<typename M1::value_type,typename M2::value_type>(),
    B,
    compose(
      plus<typename M2::value_type>(),
      bind1st(multiplies<typename M2::value_type>(),s),
      bind1st(multiplies<typename M2::value_type>(),t) 
      ),
    C
    ); 
  }

template <class M1,class M2> inline 
typename M1::value_type matrix_dot_product(const M1& A, const M2& B)
  { 
  return matrix_transform_value
    (
    0,
    incr_by<typename M1::value_type,typename M2::value_type>(),
    A,
    multiplies<typename M1::value_type>(),
    B
    ); 
  }

template <class M1,class M2,class M3> inline 
void vector_cross_product(M1& A, const M2& B, const M3& C)
  {
  invariant 
    (
    (A.rows()<=B.rows()) && (A.cols()<=B.cols()) &&
    (A.rows()<=C.rows()) && (A.cols()<=C.cols()),
    "operation assumes relaxed matrix dimensions",
    SOURCELOC
    );
  invariant 
    (
    (A.rows()==3) && (A.cols()==1),
    "operation assumes 3d vector",
    SOURCELOC
    );
  A[0] =   B[1]*C[2]-B[2]*C[1] ;
  A[1] = -(B[0]*C[2]-B[2]*C[0]);
  A[2] =   B[0]*C[1]-B[1]*C[0] ;
  }

template <class M1,class M2,class M3> inline 
void matrix_set_product(M1& R, const M2& A, const M3& B)
  {
  set_to<typename M1::value_type,float> assign_op;
  invariant( (A.cols()==B.rows()) , "matrix_set_product: dimension mismatch");
  int r,c,i;
  TData h;
  // Gr"o\3e anpassen
  R.adjust(A.rows(),B.cols());
  // Multiplikation A=B*C ausf"uhren
  for (r=0;r<R.rows();r++)
    // HOLD n=c+r*col
    for (c=0;c<R.cols();c++)
      {
      h=0;
      // i-te Spalte (row) von A mit i-ter Zeile (col) von B multiplizieren
      for (i=0;i<A.cols();i++)
        h+=(TData)(A(r,i)*B(i,c));
      assign_op(R(r,c),h);
      }
  }

template <class M1,class M2,class M3> inline 
void matrix_inc_scaled_product(M1& R,const typename M1::value_type s,const M2& A, const M3& B)
  {
  incr_by<typename M1::value_type,float> incr_op;
  invariant( (A.cols()==B.rows()) , "matrix_set_product: dimension mismatch");
  int r,c,i;
  TData h;
  // Gr"o\3e anpassen
  R.adjust(A.rows(),B.cols());
  // Multiplikation A=B*C ausf"uhren
  for (r=0;r<R.rows();r++)
    // HOLD n=c+r*col
    for (c=0;c<R.cols();c++)
      {
      h=0;
      // i-te Spalte (row) von A mit i-ter Zeile (col) von B multiplizieren
      for (i=0;i<A.cols();i++)
        h+=(TData)(A(r,i)*B(i,c));
      incr_op(R(r,c),s*h);
      }
  }
template <class M> inline
typename M::value_type matrix_det(const M& A)
  {
  typename M::value_type d(0);
  invariant(A.rows()==A.cols(),"det(nxn matrix)",SOURCELOC);

  if (A.rows()==2)
    {
    d = A(0,0)*A(1,1) - A(0,1)*A(1,0);
    }
  else
    {
    DynMatrix<0,0,typename M::value_type> buffer;
    buffer.adjust(A.rows(),A.cols());
    matrix_set(buffer,A);
    typename M::value_type factor,pivot;
    int col,row,i;

    for (col=0;col<A.cols();++col)
      {
      pivot = buffer(col,col);
      if (pivot==0) { return 0; }
      else
        {
        for (row=col+1;row<A.rows();++row)
          {
          factor = -buffer(row,col) / pivot;
          for (i=0;i<A.cols();++i)
            {
            buffer(row,i) += buffer(col,i) * factor;
            }
          }
        }
      }

    d = 1;
    for (col=0;col<A.cols();++col)
      { d *= buffer(col,col); }
    }
  return d;
  }

template <class M1> inline 
typename M1::value_type matrix_square_norm(const M1& A)
  { 
  return matrix_transform_value
    (
    0,
    incr_by<typename M1::value_type,typename M1::value_type>(),
    A,
    multiplies<typename M1::value_type>(),
    A
    ); 
  }

template <class M1,class M2> inline 
typename M1::value_type matrix_maximum_distance(const M1& A, const M2& B)
  { 
  return matrix_transform_value
    (
    0,
    set_to_max<typename M1::value_type,typename M1::value_type>(),
    A,
    difference_abs<typename M1::value_type,typename M2::value_type>(),
    B
    ); 
  }

template <class M1,class M2> inline 
typename M1::value_type matrix_square_distance(const M1& A, const M2& B)
  { 
  return matrix_transform_value
    (
    0,
    incr_by<typename M1::value_type,typename M1::value_type>(),
    A,
    difference_sqr<typename M1::value_type,typename M2::value_type>(),
    B
    ); 
  }

template <class V> inline 
void 
vector_orthogonalize
  (
  V& A
  )        
  {
  invariant
    (
    A.rows() > 1,
    "orthogonalize assumes vector (rows>1)",
    SOURCELOC
    );
  invariant
    (
    A.cols() == 1,
    "orthogonalize assumes vector (cols==1)",
    SOURCELOC
    );

  int non_zero_axis[2];
  int non_zero_counter = 0;
  int axis = 0;
  while (    (non_zero_counter<2) 
          && (axis<A.rows())
        )
    {
    if (A(axis)!=0.0)
      {
      non_zero_axis[non_zero_counter] = axis;
      ++non_zero_counter;
      }
    ++axis;
    }
  invariant
    (
    non_zero_counter!=0,
    "zero-vector cannot be orthogonalized",
    SOURCELOC
    );
  if (non_zero_counter==1)
    {
    axis = non_zero_axis[0];
    if (axis==0)
      {
      A(axis+1) = 1.0; // Assumption: A.rows()>2, see invariant above
      }
    else
      {
      A(axis-1) = 1.0; // axis is not the first axis since axis>0
      }
    A(axis) = 0.0;
    }
  else // non_zero_counter>1
    {
    swap(A(non_zero_axis[0]),A(non_zero_axis[1]));
    A(non_zero_axis[0]) *= -1.0;
    for (axis=0;axis<A.rows();++axis)
      {
      if ((axis!=non_zero_axis[0]) && (axis!=non_zero_axis[1]))
        {
        A(axis) = 0.0;
        }
      }
    }
  }

template <class M1> inline 
void matrix_normalize(M1& A)
  {
  typename M1::value_type n(sqrt(matrix_square_norm(A)));
  if (n!=0) matrix_scale(A,1.0/sqrt(matrix_square_norm(A)));
  }
template <class M> inline 
void write_matrix(ostream& os,const M& A)
  {
  os << "(";
    for (int c=0;c<A.cols();++c)
      {
      if (A.cols()>1) os << "("; 
      for (int r=0;r<A.rows();++r)
        {
        os << A(r,c);
        if (r<A.rows()-1) os << " ";
        }
      if (A.cols()>1) os << ")"; 
      }
  os << ")";
  }
template <class M> inline 
bool read_matrix(istream& is, M& A)
  {
  bool reshape(false);
  read_white(is);
  if (!is_followed_by(is,"(",false)) SYNTAX_ERROR_EXPECTED("(");
  if (is_followed_by(is,"("))
    { // matrix
    int c(0);
    int mrs(0);
    while (is_followed_by(is,"(",false))
      {
      if (c+1>A.cols()) { reshape=true; A.adjust(A.rows(),c+1); }
      int r(0);
      while (!is_followed_by(is,")"))
        {
        if (r+1>A.rows()) { reshape=true; A.adjust(r+1,A.cols()); }
        is >> A(r,c);
        ++r;
        }
      if (r>mrs) mrs=r;
      ++c;
      if (!is_followed_by(is,")",false)) SYNTAX_ERROR_EXPECTED(")");
      }
    if ((c<A.cols()) || (mrs<A.rows())) { reshape=true; A.adjust(mrs,c); }
    }
  else
    { // vector
    if (A.cols()!=1) { reshape=true; A.adjust(A.rows(),1); }
    int c(0);
    int r(0);
    while (!is_followed_by(is,")"))
      {
      if (r+1>A.rows()) { reshape=true; A.adjust(r+1,A.cols()); }
      is >> A(r,c);
      ++r;
      }
    if (r<A.rows()) { reshape=true; A.adjust(r,1); }
    }
  if (!is_followed_by(is,")",false)) SYNTAX_ERROR_EXPECTED(")");
  return !reshape;
  }
template <class M> inline 
void matrix_read_colon(istream& is, M& A)
  {
  FUNCLOG("matrix_read_colon");

  int r(-1),c(0),maxc(0);
  A.adjust(5,2); // as a default 
  do
    {
    ++r;
    if (r>=A.rows()) A.adjust(r+5,maxc+1);
    is >> A(r,c);
  trace("read",A);
    while (is_followed_by(is,",",false))
      {
      ++c;
      if (c>=A.cols()) A.adjust(r+1,c+1);
      is >> A(r,c);
      }
    maxc=max(maxc,c); c=0;
    }
  while (is_followed_by(is,":",false));
  A.adjust(r+1,maxc+1);

  trace("read",A);
  }
template <class M> inline 
void matrix_write_colon(ostream& os, const M& A)
  {
  for (int r=0;r<A.rows();++r)
    {
    for (int c=0;c<A.cols();++c)
      {
      os << A(r,c);
      if (c<A.cols()-1) os << ",";
      }
    if (r<A.rows()-1) os << ":";
    }
  }



#endif /* matvecop_HEADER */