fixablestatemodel.cpp

dysii是一款非常出色的滤波函数库

#include "FixableStateModel.hpp"

#include "boost/numeric/ublas/operation.hpp"
#include "boost/numeric/ublas/operation_sparse.hpp"

namespace aux = indii::ml::aux;
namespace ublas = boost::numeric::ublas;

using namespace indii::ml::filter;

FixableStateModel::FixableStateModel() {
  //
}

FixableStateModel::FixableStateModel(const unsigned int N) : N(N), F(0),
    fixed(N), projectState(N,N) {
  unsigned int i;
  
  fixed.clear();
  projectState.clear();
  for (i = 0; i < N; i++) {
    projectState(i,i) = 1;
  }
}

FixableStateModel::~FixableStateModel() {
  //
}

unsigned int FixableStateModel::getVariableSize() const {
  return N;
}

unsigned int FixableStateModel::getFixedSize() const {
  return F;
}

void FixableStateModel::fixVariable(const unsigned int i,
    const double value) {
  /* pre-condition */
  assert (i < N + F);

  unsigned int row, col;
  bool isFixed = true; // is variable already fixed?
  
  for (row = 0; row < projectState.size1(); row++) {
    if (projectState(row,i) == 1) {
      isFixed = false;
      break;
    }
  }

  if (!isFixed) {
    /* update projection matrix */
    /**
     * @todo Preservation in resize() is not implemented in uBLAS
     * for sparse matrices. Copy into dense matrix at present,
     * sparse matrices are used here for computational rather than
     * storage efficiency after all. Review in future if
     * preservation for sparse matrices is implemented in uBLAS.
     */
    aux::matrix projectStateDense(projectState);
    ublas::range cols(0, N + F);
    ublas::range to(row, N - 1);
    ublas::range from(row + 1, N);
    ublas::matrix_range<aux::matrix>(projectStateDense, to, cols) =
        ublas::matrix_range<aux::matrix>(projectStateDense, from, cols);
    projectStateDense.resize(N - 1, N + F, true);
    projectState.resize(N - 1, N + F, false);
    projectState.clear();
    for (col = 0; col < projectStateDense.size2(); col++) {
      for (row = 0; row < projectStateDense.size1(); row++) {
    	  if (projectStateDense(row,col) == 1) {
          projectState(row,col) = 1;
    	  }
    	}
    }
  }
  
  /* update fixed variables */
  fixed(i) = value;
  N--;
  F++;

  /* post-condition */
  assert (projectState.size1() == N);
  assert (projectState.size2() == N + F);
}

aux::vector FixableStateModel::toState(const aux::vector& x) const {
  /* pre-condition */
  assert (x.size() == N + F);

  if (F == 0) {
    return x;
  } else {
    aux::vector result(N);
    ublas::axpy_prod(projectState, x, result, true);
    return result;
  }
}

aux::vector FixableStateModel::fromState(const aux::vector& x) const {
  /* pre-condition */
  assert (x.size() == N);

  if (F == 0) {
    return x;
  } else {
    aux::vector result(fixed);
    ublas::axpy_prod(x, projectState, result, false);
    return result;
  }
}

aux::symmetric_matrix FixableStateModel::toState(
    const aux::symmetric_matrix& x) const {
  /* pre-condition */
  assert (x.size1() == N + F);

  if (F == 0) {
    return x;
  } else {
    aux::symmetric_matrix result(N);
    aux::matrix tmp(N, N + F);
    ublas::axpy_prod(projectState, x, tmp, true);
    ublas::axpy_prod(tmp, trans(projectState), result, true);
    return result;
  }
}

aux::symmetric_matrix FixableStateModel::fromState(
    const aux::symmetric_matrix& x) const {
  /* pre-condition */
  assert (x.size1() == N);

  if (F == 0) {
    return x;
  } else {
    aux::symmetric_matrix result(N + F);
    aux::matrix tmp(N + F, N);
    ublas::axpy_prod(trans(projectState), x, tmp, true);
    ublas::axpy_prod(tmp, projectState, result, true);
    return result;
  }
}