ConstantStepOdeSolver Class Reference

#include <ConstantStepOdeSolver.h>

Inheritance diagram for ConstantStepOdeSolver:

Public Member Functions
	ConstantStepOdeSolver ()
	ConstantStepOdeSolver (const AbsODERHS &rhs)
	ConstantStepOdeSolver (const ConstantStepOdeSolver &r)
double	getCoordinate (const unsigned which, const unsigned idx) const
double	getIntegrated (unsigned which, unsigned idx) const
double	getPeakTime (unsigned which) const
const AbsODERHS *	getRHS () const
AbsODERHS *	getRHS ()
double	getTime (const unsigned idx) const
double	interpolateCoordinate (unsigned which, double t, bool cubic=false) const
double	interpolateIntegrated (unsigned which, double t, bool cubic=false) const
double	lastDeltaT () const
unsigned	lastDim () const
double	lastMaxT () const
unsigned	lastRunLength () const
virtual const char *	methodName () const =0
ConstantStepOdeSolver &	operator= (const ConstantStepOdeSolver &r)
void	run (const double *initialConditions, unsigned lenConditions, double dt, unsigned nSteps)
void	setHistory (double dt, const double *data, unsigned dim, unsigned runLen)
void	setRHS (const AbsODERHS &rhs)
void	truncateCoordinate (unsigned which, double minValue, double maxValue)
void	writeHistory (std::ostream &os, double dt, bool cubic=false) const
void	writeIntegrated (std::ostream &os, unsigned which, double dt, bool cubic=false) const
virtual	~ConstantStepOdeSolver ()

Protected Attributes
AbsODERHS *	rhs_

Private Member Functions
void	integrateCoordinate (const unsigned which)
virtual void	step (double t, double dt, const double *x, unsigned lenX, double *coordIncrement) const =0

Private Attributes
std::vector< double >	chargeBuffer_
unsigned	dim_
double	dt_
std::vector< double >	historyBuffer_
unsigned	lastIntegrated_
unsigned	runLen_

Detailed Description

Definition at line 14 of file ConstantStepOdeSolver.h.

Constructor & Destructor Documentation

ConstantStepOdeSolver::ConstantStepOdeSolver ( )

inline

Definition at line 17 of file ConstantStepOdeSolver.h.

        : rhs_(0), dt_(0.0), dim_(0), runLen_(0), lastIntegrated_(0) {}

ConstantStepOdeSolver::ConstantStepOdeSolver ( const AbsODERHS &  rhs )

inline

Definition at line 20 of file ConstantStepOdeSolver.h.

References AbsODERHS::clone(), and rhs_.

                                                       :
        rhs_(0), dt_(0.0), dim_(0), runLen_(0), lastIntegrated_(0)
    {
        rhs_ = rhs.clone();
    }

ConstantStepOdeSolver::ConstantStepOdeSolver

(

const ConstantStepOdeSolver &

r

)

Definition at line 145 of file ConstantStepOdeSolver.cc.

References AbsODERHS::clone(), and rhs_.

    : rhs_(0),
      dt_(r.dt_),
      dim_(r.dim_),
      runLen_(r.runLen_),
      lastIntegrated_(r.lastIntegrated_),
      historyBuffer_(r.historyBuffer_),
      chargeBuffer_(r.chargeBuffer_)
{
    if (r.rhs_)
        rhs_ = r.rhs_->clone();
}

virtual ConstantStepOdeSolver::~ConstantStepOdeSolver ( )

inlinevirtual

Definition at line 30 of file ConstantStepOdeSolver.h.

References rhs_.

{delete rhs_;}

Member Function Documentation

double ConstantStepOdeSolver::getCoordinate ( const unsigned  which, const unsigned  idx  ) const

inline

Definition at line 54 of file ConstantStepOdeSolver.h.

References dim_, edm::hlt::Exception, historyBuffer_, and runLen_.

Referenced by writeHistory().

    {
        if (which >= dim_ || idx >= runLen_) throw cms::Exception(
            "In ConstantStepOdeSolver::getCoordinate: index out of range");
        return historyBuffer_[dim_*idx + which];
    }

double ConstantStepOdeSolver::getIntegrated	(	unsigned	which,
		unsigned	idx
	)		const

Definition at line 39 of file ConstantStepOdeSolver.cc.

References chargeBuffer_, dim_, edm::hlt::Exception, customizeTrackingMonitorSeedNumber::idx, lastIntegrated_, and runLen_.

Referenced by writeIntegrated().

{
    if (which >= dim_ || idx >= runLen_) throw cms::Exception(
        "In ConstantStepOdeSolver::getIntegrated: index out of range");
    if (lastIntegrated_ != which)
        (const_cast<ConstantStepOdeSolver*>(this))->integrateCoordinate(which);
    return chargeBuffer_[idx];
}

double ConstantStepOdeSolver::getPeakTime

(

unsigned

which

)

const

Definition at line 11 of file ConstantStepOdeSolver.cc.

References dim_, dt_, edm::hlt::Exception, historyBuffer_, i, prof2calltree::l, alignCSCRings::r, and runLen_.

Referenced by QIE8Simulator::preampPeakTime().

{
    if (which >= dim_) throw cms::Exception(
        "In ConstantStepOdeSolver::getPeakTime: index out of range");
    if (runLen_ < 3) throw cms::Exception(
        "In ConstantStepOdeSolver::getPeakTime: not enough data");

    const double* hbuf = &historyBuffer_[which];
    double maxval = hbuf[0];
    unsigned maxind = 0;
    for (unsigned i=1; i<runLen_; ++i)
        if (hbuf[dim_*i] > maxval)
        {
            maxval = hbuf[dim_*i];
            maxind = i;
        }
    if (maxind == 0U)
        return 0.0;
    if (maxind == runLen_-1U)
        return dt_*maxind;
    const double l = hbuf[dim_*(maxind - 1U)];
    const double r = hbuf[dim_*(maxind + 1U)];
    if (l < maxval || r < maxval)
        return dt_*(maxind + (l - r)/2.0/(l + r - 2.0*maxval));
    else
        return dt_*maxind;
}

const AbsODERHS* ConstantStepOdeSolver::getRHS ( ) const

inline

Definition at line 38 of file ConstantStepOdeSolver.h.

References rhs_.

Referenced by QIE8Simulator::getRHS(), and QIE8Simulator::modifiableRHS().

{return rhs_;}

AbsODERHS* ConstantStepOdeSolver::getRHS ( )

inline

Definition at line 39 of file ConstantStepOdeSolver.h.

References rhs_.

{return rhs_;}

double ConstantStepOdeSolver::getTime ( const unsigned  idx ) const

inline

Definition at line 47 of file ConstantStepOdeSolver.h.

References dt_, edm::hlt::Exception, and runLen_.

Referenced by writeHistory(), and writeIntegrated().

    {
        if (idx >= runLen_) throw cms::Exception(
            "In ConstantStepOdeSolver::getTime: index out of range");
        return idx*dt_;
    }

void ConstantStepOdeSolver::integrateCoordinate ( const unsigned  which )

private

Definition at line 49 of file ConstantStepOdeSolver.cc.

References chargeBuffer_, dim_, dt_, edm::hlt::Exception, historyBuffer_, i, lastIntegrated_, and runLen_.

{
    if (runLen_ < 4) throw cms::Exception(
        "In ConstantStepOdeSolver::integrateCoordinate: not enough data");
    if (chargeBuffer_.size() < runLen_)
        chargeBuffer_.resize(runLen_);
    double* integ = &chargeBuffer_[0];
    const double* coord = &historyBuffer_[which];

    integ[0] = 0.0;
    integ[1] = coord[dim_*0]*(3.0/8.0) +
               coord[dim_*1]*(19.0/24.0) +
               coord[dim_*2]*(-5.0/24.0) +
               coord[dim_*3]*(1.0/24.0);
    long double sum = integ[1];
    const unsigned rlenm1 = runLen_ - 1U;
    for (unsigned i=2; i<rlenm1; ++i)
    {
        sum += (coord[dim_*(i-2U)]*(-1.0/24.0) +
                coord[dim_*(i-1U)]*(13.0/24.0) +
                coord[dim_*i]     *(13.0/24.0) +
                coord[dim_*(i+1U)]*(-1.0/24.0));
        integ[i] = sum;
    }
    sum += (coord[dim_*rlenm1]     *(3.0/8.0) +
            coord[dim_*(rlenm1-1U)]*(19.0/24.0) +
            coord[dim_*(rlenm1-2U)]*(-5.0/24.0) +
            coord[dim_*(rlenm1-3U)]*(1.0/24.0));
    integ[rlenm1] = sum;                                      
    if (dt_ != 1.0)
        for (unsigned i=1; i<runLen_; ++i)
            integ[i] *= dt_;
    lastIntegrated_ = which;
}

double ConstantStepOdeSolver::interpolateCoordinate	(	unsigned	which,
		double	t,
		bool	cubic = false
	)		const

Definition at line 195 of file ConstantStepOdeSolver.cc.

References newFWLiteAna::base, dim_, dt_, edm::hlt::Exception, historyBuffer_, interpolateLinear(), runLen_, and x.

Referenced by QIE8Simulator::controlOutput(), QIE8Simulator::getCharge(), QIE8Simulator::preampOutput(), and writeHistory().

{
    if (which >= dim_) throw cms::Exception(
        "In ConstantStepOdeSolver::interpolateCoordinate: index out of range");
    if (runLen_ < 2U || (cubic && runLen_ < 4U)) throw cms::Exception(
        "In ConstantStepOdeSolver::interpolateCoordinate: not enough data");
    const double maxt = runLen_ ? dt_*(runLen_-1U) : 0.0;
    if (t < 0.0 || t > maxt) throw cms::Exception(
        "In ConstantStepOdeSolver::interpolateCoordinate: time out of range");

    const double* arr = &historyBuffer_[0];
    if (t == 0.0)
        return arr[which];
    else if (t == maxt)
        return arr[which + dim_*(runLen_-1U)];

    // Translate time into timestep units
    const double tSteps = t/dt_;
    unsigned nLow = tSteps;
    if (nLow >= runLen_ - 1)
        nLow = runLen_ - 2;
    double x = tSteps - nLow;

    if (cubic)
    {
        unsigned i0 = 0;
        if (nLow == runLen_ - 2)
        {
            i0 = nLow - 2U;
            x += 2.0;
        }
        else if (nLow)
        {
            i0 = nLow - 1U;
            x += 1.0;
        }
        const double* base = arr + (which + dim_*i0);
        return interpolateLinear(x*(3.0 - x)/2.0,
                                 interpolateLinear(x/3.0, base[0], base[dim_*3]),
                                 interpolateLinear(x-1.0, base[dim_], base[dim_*2]));
    }
    else
        return interpolateLinear(x, arr[which+dim_*nLow], arr[which+dim_*(nLow+1U)]);
}

double ConstantStepOdeSolver::interpolateIntegrated	(	unsigned	which,
		double	t,
		bool	cubic = false
	)		const

Definition at line 241 of file ConstantStepOdeSolver.cc.

References newFWLiteAna::base, chargeBuffer_, dim_, dt_, edm::hlt::Exception, interpolateLinear(), lastIntegrated_, runLen_, and x.

Referenced by QIE8Simulator::getCharge(), and writeIntegrated().

{
    if (which >= dim_) throw cms::Exception(
        "In ConstantStepOdeSolver::interpolateIntegrated: index out of range");
    if (runLen_ < 2U || (cubic && runLen_ < 4U)) throw cms::Exception(
        "In ConstantStepOdeSolver::interpolateIntegrated: not enough data");
    const double maxt = runLen_ ? dt_*(runLen_-1U) : 0.0;
    if (t < 0.0 || t > maxt) throw cms::Exception(
        "In ConstantStepOdeSolver::interpolateIntegrated: time out of range");
    if (lastIntegrated_ != which)
        (const_cast<ConstantStepOdeSolver*>(this))->integrateCoordinate(which);

    const double* buf = &chargeBuffer_[0];
    if (t == 0.0)
        return buf[0];
    else if (t == maxt)
        return buf[runLen_-1U];

    // Translate time into timestep units
    const double tSteps = t/dt_;
    unsigned nLow = tSteps;
    if (nLow >= runLen_ - 1)
        nLow = runLen_ - 2;
    double x = tSteps - nLow;

    if (cubic)
    {
        unsigned i0 = 0;
        if (nLow == runLen_ - 2)
        {
            i0 = nLow - 2U;
            x += 2.0;
        }
        else if (nLow)
        {
            i0 = nLow - 1U;
            x += 1.0;
        }
        const double* base = buf + i0;
        return interpolateLinear(x*(3.0 - x)/2.0,
                                 interpolateLinear(x/3.0, base[0], base[3]),
                                 interpolateLinear(x-1.0, base[1], base[2]));
    }
    else
        return interpolateLinear(x, buf[nLow], buf[nLow+1U]);
}

double ConstantStepOdeSolver::lastDeltaT ( ) const

inline

Definition at line 44 of file ConstantStepOdeSolver.h.

References dt_.

Referenced by operator<<().

{return dt_;}

unsigned ConstantStepOdeSolver::lastDim ( ) const

inline

Definition at line 42 of file ConstantStepOdeSolver.h.

References dim_.

{return dim_;}

double ConstantStepOdeSolver::lastMaxT ( ) const

inline

Definition at line 45 of file ConstantStepOdeSolver.h.

References dt_, and runLen_.

Referenced by QIE8Simulator::lastStopTime(), QIE8Simulator::run(), writeHistory(), and writeIntegrated().

{return runLen_ ? dt_*(runLen_-1U) : 0.0;}

unsigned ConstantStepOdeSolver::lastRunLength ( ) const

inline

Definition at line 43 of file ConstantStepOdeSolver.h.

References runLen_.

{return runLen_;}

virtual const char* ConstantStepOdeSolver::methodName ( ) const

pure virtual

Implemented in RK4, RK2, and EulerOdeSolver.

Referenced by writeHistory(), and writeIntegrated().

ConstantStepOdeSolver & ConstantStepOdeSolver::operator=

(

const ConstantStepOdeSolver &

r

)

Definition at line 158 of file ConstantStepOdeSolver.cc.

References chargeBuffer_, AbsODERHS::clone(), dim_, dt_, historyBuffer_, lastIntegrated_, rhs_, and runLen_.

{
    if (this != &r)
    {
        delete rhs_; rhs_ = 0;
        dt_ = r.dt_;
        dim_ = r.dim_;
        runLen_ = r.runLen_;
        lastIntegrated_ = r.lastIntegrated_;
        historyBuffer_ = r.historyBuffer_;
        chargeBuffer_ = r.chargeBuffer_;
        if (r.rhs_)
            rhs_ = r.rhs_->clone();
    }
    return *this;
}

void ConstantStepOdeSolver::run	(	const double *	initialConditions,
		unsigned	lenConditions,
		double	dt,
		unsigned	nSteps
	)

Definition at line 312 of file ConstantStepOdeSolver.cc.

References assert(), dim_, dt, dt_, edm::hlt::Exception, historyBuffer_, i, lastIntegrated_, GetRecoTauVFromDQM_MC_cff::next, rhs_, runLen_, step(), and tree::t.

Referenced by Types.EventID::cppID(), Types.LuminosityBlockID::cppID(), and QIE8Simulator::run().

{
    if (!nSteps || !lenInitialConditions)
        return;
    if (!rhs_) throw cms::Exception(
        "In ConstantStepOdeSolver::run: ODE right hand side has not been set");
    if (dt <= 0.0) throw cms::Exception(
        "In ConstantStepOdeSolver::run: can not run backwards in time");
    assert(initialConditions);
    dt_ = dt;
    dim_ = lenInitialConditions;
    lastIntegrated_ = dim_;
    runLen_ = nSteps + 1U;
    const unsigned sz = (runLen_+1U)*dim_;
    if (historyBuffer_.size() < sz)
        historyBuffer_.resize(sz);
    double* arr = &historyBuffer_[0];
    for (unsigned i=0; i<lenInitialConditions; ++i)
        arr[i] = initialConditions[i];
    double* stepBuffer = arr + runLen_*dim_;

    for (unsigned i = 0; i < nSteps; ++i, arr += lenInitialConditions)
    {
        const double t = i*dt;
        this->step(t, dt, arr, lenInitialConditions, stepBuffer);
        double* next = arr + lenInitialConditions;
        for (unsigned i=0; i<lenInitialConditions; ++i)
            next[i] = arr[i] + stepBuffer[i];
    }
}

void ConstantStepOdeSolver::setHistory	(	double	dt,
		const double *	data,
		unsigned	dim,
		unsigned	runLen
	)

Definition at line 289 of file ConstantStepOdeSolver.cc.

References assert(), dim_, dt, dt_, edm::hlt::Exception, historyBuffer_, i, lastIntegrated_, and runLen_.

Referenced by QIE8Simulator::run().

{
    const unsigned len = dim*runLen;
    if (!len)
        return;
    if (dt <= 0.0) throw cms::Exception(
        "In ConstantStepOdeSolver::setHistory: can not run backwards in time");
    assert(data);
    const unsigned sz = dim*(runLen + 1U);
    if (historyBuffer_.size() < sz)
        historyBuffer_.resize(sz);
    dt_ = dt;
    dim_ = dim;
    runLen_ = runLen;
    lastIntegrated_ = dim_;
    double* arr = &historyBuffer_[0];
    for (unsigned i=0; i<len; ++i)
        *arr++ = *data++;
    for (unsigned i=0; i<dim; ++i)
        *arr++ = 0.0;
}

void ConstantStepOdeSolver::setRHS ( const AbsODERHS &  rhs )

inline

Definition at line 33 of file ConstantStepOdeSolver.h.

References AbsODERHS::clone(), and rhs_.

Referenced by QIE8Simulator::setRHS().

    {
        delete rhs_;
        rhs_ = rhs.clone();
    }

virtual void ConstantStepOdeSolver::step ( double  t, double  dt, const double *  x, unsigned  lenX, double *  coordIncrement  ) const

privatepure virtual

Implemented in RK4, RK2, and EulerOdeSolver.

Referenced by run().

void ConstantStepOdeSolver::truncateCoordinate	(	unsigned	which,
		double	minValue,
		double	maxValue
	)

Definition at line 176 of file ConstantStepOdeSolver.cc.

References dim_, edm::hlt::Exception, historyBuffer_, i, and runLen_.

Referenced by QIE8Simulator::run().

{
    if (which >= dim_) throw cms::Exception(
        "In ConstantStepOdeSolver::truncateCoordinate: index out of range");
    if (minValue > maxValue) throw cms::Exception(
        "In ConstantStepOdeSolver::truncateCoordinate: invalid truncation range");

    double* buf = &historyBuffer_[which];
    for (unsigned i=0; i<runLen_; ++i)
    {
        if (buf[dim_*i] < minValue)
            buf[dim_*i] = minValue;
        else if (buf[dim_*i] > maxValue)
            buf[dim_*i] = maxValue;
    }
}

void ConstantStepOdeSolver::writeHistory	(	std::ostream &	os,
		double	dt,
		bool	cubic = false
	)		const

Definition at line 84 of file ConstantStepOdeSolver.cc.

References dim_, dt, dt_, getCoordinate(), getTime(), i, interpolateCoordinate(), lastMaxT(), methodName(), runLen_, tree::t, and tmax.

Referenced by operator<<().

{
    os << "# " << this->methodName() << '\n';
    if (dim_ && runLen_)
    {
        if (dt == dt_)
        {
            for (unsigned ipt=0; ipt<runLen_; ++ipt)
            {
                os << getTime(ipt);
                for (unsigned which=0; which<dim_; ++which)
                    os << ' ' << getCoordinate(which, ipt);
                os << '\n';
            }
        }
        else
        {
            const double tmax = lastMaxT();
            for (unsigned i=0; ; ++i)
            {
                const double t = i*dt;
                if (t > tmax)
                    break;
                os << t;
                for (unsigned which=0; which<dim_; ++which)
                    os << ' ' << interpolateCoordinate(which, t, cubic);
                os << '\n';
            }
        }
    }
}

void ConstantStepOdeSolver::writeIntegrated	(	std::ostream &	os,
		unsigned	which,
		double	dt,
		bool	cubic = false
	)		const

Definition at line 118 of file ConstantStepOdeSolver.cc.

References dim_, dt, dt_, getIntegrated(), getTime(), i, interpolateIntegrated(), lastMaxT(), methodName(), runLen_, tree::t, and tmax.

{
    os << "# " << this->methodName() << '\n';
    if (dim_ && runLen_)
    {
        if (dt == dt_)
        {
            for (unsigned ipt=0; ipt<runLen_; ++ipt)
                os << getTime(ipt) << ' ' << getIntegrated(which, ipt) << '\n';
        }
        else
        {
            const double tmax = lastMaxT();
            for (unsigned i=0; ; ++i)
            {
                const double t = i*dt;
                if (t > tmax)
                    break;
                os << t << ' ' << interpolateIntegrated(which, t, cubic) << '\n';
            }
        }
    }
}

Member Data Documentation

std::vector<double> ConstantStepOdeSolver::chargeBuffer_

private

Definition at line 118 of file ConstantStepOdeSolver.h.

Referenced by getIntegrated(), integrateCoordinate(), interpolateIntegrated(), and operator=().

unsigned ConstantStepOdeSolver::dim_

private

Definition at line 113 of file ConstantStepOdeSolver.h.

Referenced by getCoordinate(), getIntegrated(), getPeakTime(), integrateCoordinate(), interpolateCoordinate(), interpolateIntegrated(), lastDim(), operator=(), run(), setHistory(), truncateCoordinate(), writeHistory(), and writeIntegrated().

double ConstantStepOdeSolver::dt_

private

Definition at line 112 of file ConstantStepOdeSolver.h.

Referenced by getPeakTime(), getTime(), integrateCoordinate(), interpolateCoordinate(), interpolateIntegrated(), lastDeltaT(), lastMaxT(), operator=(), run(), setHistory(), writeHistory(), and writeIntegrated().

std::vector<double> ConstantStepOdeSolver::historyBuffer_

private

Definition at line 117 of file ConstantStepOdeSolver.h.

Referenced by getCoordinate(), getPeakTime(), integrateCoordinate(), interpolateCoordinate(), operator=(), run(), setHistory(), and truncateCoordinate().

unsigned ConstantStepOdeSolver::lastIntegrated_

private

Definition at line 115 of file ConstantStepOdeSolver.h.

Referenced by getIntegrated(), integrateCoordinate(), interpolateIntegrated(), operator=(), run(), and setHistory().

AbsODERHS* ConstantStepOdeSolver::rhs_

protected

Definition at line 100 of file ConstantStepOdeSolver.h.

Referenced by ConstantStepOdeSolver(), getRHS(), operator=(), run(), setRHS(), EulerOdeSolver::step(), RK2::step(), RK4::step(), and ~ConstantStepOdeSolver().

unsigned ConstantStepOdeSolver::runLen_

private

Definition at line 114 of file ConstantStepOdeSolver.h.

Referenced by getCoordinate(), getIntegrated(), getPeakTime(), getTime(), integrateCoordinate(), interpolateCoordinate(), interpolateIntegrated(), lastMaxT(), lastRunLength(), operator=(), run(), setHistory(), truncateCoordinate(), writeHistory(), and writeIntegrated().