Logo Search packages:      
Sourcecode: kboincspy version File versions  Download package

kbslhcinterpolator.cpp

/***************************************************************************
 *   Copyright (C) 2005 by Roberto Virga                                   *
 *   rvirga@users.sf.net                                                   *
 *                                                                         *
 *   This program 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; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program 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, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include <math.h>

#include <kbslhctaskmonitor.h>

#include "kbslhcinterpolator.h"

KBSLHCInterpolator::KBSLHCInterpolator(unsigned set, unsigned position,
                                       KBSLHCTaskMonitor *parent, const char *name)
                  : QObject(parent, name), m_set(set), m_position(position)
{
  resetIndices();
  update();
  
  connect(parent, SIGNAL(updatedSet(unsigned)), this, SLOT(update(unsigned)));  
}

unsigned KBSLHCInterpolator::set() const
{
  return m_set;
}

unsigned KBSLHCInterpolator::position() const
{
  return m_position;
}

double KBSLHCInterpolator::interpolateXCoord(double turn)
{
  if(m_indices.isEmpty())
    return 0.0;
  else if(turn <= m_indices.first())
    return m_data[m_indices.first()].coord.x.mm;
  else if(turn >= m_indices.last())
    return m_data[m_indices.last()].coord.x.mm;
    
  computeIndices(turn);
  computeCoefficients(turn);
  
  double output = 0.0;
  for(int i = m_start; i <= m_end; ++i)
    output += m_data[*m_index[i]].coord.x.mm * m_numerator[i] / m_denominator[i];
  
  return output;
}

double KBSLHCInterpolator::interpolateYCoord(double turn)
{
  if(m_indices.isEmpty())
    return 0.0;
  else if(turn <= m_indices.first())
    return m_data[m_indices.first()].coord.y.mm;
  else if(turn >= m_indices.last())
    return m_data[m_indices.last()].coord.y.mm;
    
  computeIndices(turn);
  computeCoefficients(turn);
  
  double output = 0.0;
  for(int i = m_start; i <= m_end; ++i)
    output += m_data[*m_index[i]].coord.y.mm * m_numerator[i] / m_denominator[i];
  
  return output;
}

double KBSLHCInterpolator::interpolateEnergy(double turn)
{
  if(m_indices.isEmpty())
    return 0.0;
  else if(turn <= m_indices.first())
    return m_data[m_indices.first()].coord.egy;
  else if(turn >= m_indices.last())
    return m_data[m_indices.last()].coord.egy;
    
  computeIndices(turn);
  computeCoefficients(turn);
  
  double output = 0.0;
  for(int i = m_start; i <= m_end; ++i)
    output += m_data[*m_index[i]].coord.egy * m_numerator[i] / m_denominator[i];
  
  return output;
}

KBSLHCTaskMonitor *KBSLHCInterpolator::taskMonitor()
{
  return static_cast<KBSLHCTaskMonitor*>(parent());
}

void KBSLHCInterpolator::resetIndices()
{
  m_index[PrevInf] = m_index[Inf] = m_index[Sup] = m_index[SuccSup] = m_indices.constBegin();
  m_start = m_end = -1;
}

void KBSLHCInterpolator::computeIndices(double turn)
{
  if(m_indices.isEmpty()) {
    resetIndices();
    return;
  }
  
  bool updated = false;
  
  if(m_index[Sup] != m_indices.constEnd() && turn >= *m_index[Sup])
  {
    updated = true;
    do ++m_index[Sup];
    while(m_index[Sup] != m_indices.constBegin() && turn >= *m_index[Sup]);
    m_index[Inf] = m_index[Sup];
    if(m_index[Inf] != m_indices.constBegin()) --m_index[Inf];
  } else if(m_indices.constBegin() != m_index[Inf] && turn < *m_index[Inf]) {
    updated = true;
    do --m_index[Inf];
    while(m_index[Inf] != m_indices.constBegin() && turn < *m_index[Inf]);
    m_index[Sup] = m_index[Inf];
    ++m_index[Sup];
  }
  
  QValueList<unsigned>::const_iterator index;
  
  index = m_index[Inf];
  if(index != m_indices.constBegin()) --index;
  if(updated || index != m_index[PrevInf]) {
    updated = true;
    m_index[PrevInf] = index;
  }
  
  index = m_index[Sup];
  if(index != m_indices.constEnd()) ++index;
  if(updated || index != m_index[SuccSup]) {
    updated = true;
    m_index[SuccSup] = index;
  }
  
  if(updated) m_start = m_end = -1;
}
  
void KBSLHCInterpolator::computeCoefficients(double turn)
{
  if(m_start < 0)
  {
    m_end = MaxDegree-1;
    while(m_end >= 0 && m_indices.constEnd() == m_index[m_end])
      m_end--;
    
    m_start = 0;
    while(m_start < m_end && m_index[m_start] == m_index[m_start+1])
      m_start++;
    
    if(m_end >= 0)
    {
      double m_binomial[m_end+1][m_end+1];
      
      for(int i = m_start+1; i <= m_end; ++i)
        for(int j = m_start; j < i; ++j)
          m_binomial[i][j] = *m_index[i] - *m_index[j];
          
      for(int i = m_start; i <= m_end; ++i) {
        m_denominator[i] = 1.0;
        for(int j = m_start; j <= m_end; ++j)
          if(i != j) m_denominator[i] *= (i > j) ? m_binomial[i][j] : -m_binomial[j][i];
      }
    }
    
    for(int i = 0; i < m_start; ++i)
      m_denominator[i] = 0.0;
    for(int i = m_end+1; i < MaxDegree; ++i)
      m_denominator[i] = 0.0;
     
    m_turn = -1.0;
  }
  
  if(fabs(m_turn - turn) >= 1e-3)
  {
    if(m_end >= 0)
    {
      double m_binomial[m_end+1];
      
      for(int i = m_start; i <= m_end; ++i)
        m_binomial[i] = turn - *m_index[i];
        
      for(int i = m_start; i <= m_end; ++i) {
        m_numerator[i] = 1.0;
        for(int j = m_start; j <= m_end; ++j)
          if(i != j) m_numerator[i] *= m_binomial[j];
      }
    }
    
    for(int i = 0; i < m_start; ++i)
      m_numerator[i] = 0.0;
    for(int i = m_end+1; i < MaxDegree; ++i)
      m_numerator[i] = 0.0;
      
    m_turn = turn;
  }
}

void KBSLHCInterpolator::update()
{
  double query = -1.0;
  if(m_index[Inf] != m_indices.constEnd())
    query = *m_index[Inf];
  
  m_data.clear();
  m_indices.clear();
  resetIndices();
  
  const LHCState *state = taskMonitor()->state();
  if(NULL == state) return;
  
  if(!state->output.contains(m_set)) return;
  const LHCOutput *output = &state->output[m_set];
  
  const LHCHeader *header = &output->header;
  
  const unsigned count = header->part.last - header->part.fst + 1;
  if(m_position >= count) return;
  
  m_data = output->data[m_position];
  m_indices = m_data.keys();
  qHeapSort(m_indices);
  resetIndices();
  
  if(query >= 0.0)
    computeIndices(query);  
}

void KBSLHCInterpolator::update(unsigned set)
{
  if(this->set() != set) return;
  
  update();
}

#include "kbslhcinterpolator.moc"

Generated by  Doxygen 1.6.0   Back to index