hessian.hpp

/*
  hessian matrix construction code
 
*/
 
 
/*
  This file is part of LOOS.
 
  LOOS (Lightweight Object-Oriented Structure library)
  Copyright (c) 2009 Tod D. Romo
  Department of Biochemistry and Biophysics
  School of Medicine & Dentistry, University of Rochester
 
  This package (LOOS) is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation under version 3 of the License.
 
  This package is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
*/
 
 
#if !defined(LOOS_HESSIAN_HPP)
#define LOOS_HESSIAN_HPP
 
 
#include <loos.hpp>
 
#include "spring_functions.hpp"
 
namespace ENM {
 
 
  class SuperBlock {
  public:
    SuperBlock() : springs(0) { }
  
 
    SuperBlock(SpringFunction* func, const loos::AtomicGroup& nodelist) : springs(func), nodes(nodelist) { }
    SuperBlock(const SuperBlock& b) : springs(b.springs), nodes(b.nodes) { }
    virtual ~SuperBlock() { }
 
    uint size() const { return(static_cast<uint>(nodes.size())); }
 
    // ------------------------------------------------------
    virtual SpringFunction::Params setParams(const SpringFunction::Params& v) {
      return(springs->setParams(v));
    }
 
    virtual bool validParams() const { return(springs->validParams()); }
 
    virtual uint paramSize() const { return(springs->paramSize()); }
    // ------------------------------------------------------
 
    virtual loos::DoubleMatrix block(const uint j, const uint i) {
      return(blockImpl(j, i, springs));
    }
 
 
  protected:
 
 
    loos::DoubleMatrix blockImpl(const uint j, const uint i, SpringFunction* fptr) {
      if (i >= size() || j >= size())
        throw(std::runtime_error("Invalid index in Hessian SuperBlock"));
 
      if (fptr == 0)
        throw(std::runtime_error("No spring function defined for hessian!"));
 
      loos::GCoord u = nodes[j]->coords();
      loos::GCoord v = nodes[i]->coords();
      loos::GCoord d = v - u;
    
      loos::DoubleMatrix K = fptr->constant(u, v, d);
      loos::DoubleMatrix B(3, 3);
      for (uint y=0; y<3; ++y)
        for (uint x=0; x<3; ++x)
          B(x, y) = d[x]*d[y] * K(x,y);
 
      return(B);
    }
 
 
    SpringFunction* springs;
    loos::AtomicGroup nodes;
  };
 
 
 
 
  class SuperBlockDecorator : public SuperBlock {
  public:
 
    SuperBlockDecorator(SuperBlock* b) : SuperBlock(*b), decorated(b) { }
 
  protected:
    SuperBlock *decorated;
  };
 
 
 
 
  class BoundSuperBlock : public SuperBlockDecorator {
  public:
    BoundSuperBlock(SuperBlock* b, SpringFunction* bs, loos::Math::Matrix<int>& cm) :
      SuperBlockDecorator(b),
      bound_spring(bs),
      connectivity(cm)
    {
      if (connectivity.rows() != connectivity.cols() && connectivity.cols() != size())
        throw(std::runtime_error("Connectivity matrix and Nodelist have differing sizes"));
    }
 
 
    // Block now checks to see if nodes i and j are connected and, if
    // so, uses our alternative spring function.
    loos::DoubleMatrix block(const uint j, const uint i) {
      if (connectivity(j, i))
        return(blockImpl(j, i, bound_spring));
      else
        return(decorated->block(j, i));
    }
 
    SpringFunction::Params setParams(const SpringFunction::Params& v) {
      SpringFunction::Params u = bound_spring->setParams(v);
      if (! u.empty())
        u = decorated->setParams(u);
      return(u);
    }
 
    bool validParams() const { return(bound_spring->validParams() && decorated->validParams()); }
 
    uint paramSize() const { return(bound_spring->paramSize() + decorated->paramSize()); }
 
  private:
    SpringFunction* bound_spring;
    loos::Math::Matrix<int> connectivity;
  };
 
};
 
#endif