00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #ifndef __RVL_TOOLS_POLYOP_H_
00034 #define __RVL_TOOLS_POLYOP_H_
00035
00036 #include "functions.hh"
00037 #include "op.hh"
00038
00039 namespace RVL {
00040
00047 template<class Scalar>
00048 class PolynomialOperator : public OperatorWithInvertibleDeriv<Scalar> {
00049
00050 protected:
00051 Space<Scalar> & spc;
00052 std::valarray<Scalar> coef;
00053
00059 virtual void apply(const Vector<Scalar> & x,
00060 Vector<Scalar> & y) const {
00061 Vector<Scalar> temp(x);
00062 ElementwiseMultiply<Scalar> dotstar;
00063 int size = coef.size();
00064 if( size > 0) {
00065 RVLAssignConst<Scalar> constantterm(coef[0]);
00066 y.eval(constantterm);
00067 }
00068 if( size > 1 )
00069 y.linComb(coef[1], temp);
00070 for( int i = 2; i < size; i++) {
00071 temp.eval(dotstar, temp, x);
00072 y.linComb(coef[i], temp);
00073 }
00074 }
00075
00077 virtual void applyDeriv(const Vector<Scalar> & x,
00078 const Vector<Scalar> & dx,
00079 Vector<Scalar> & dy) const {
00080 Vector<Scalar> temp(x);
00081 ElementwiseMultiply<Scalar> dotstar;
00082 int size = coef.size();
00083 if( size > 1) {
00084 RVLAssignConst<Scalar> constantterm(coef[1]);
00085 dy.eval(constantterm);
00086 } else {
00087 dy.zero();
00088 }
00089 if( size > 2 )
00090 dy.linComb(coef[2]*Scalar(2), temp);
00091 for( int i = 3; i < size; i++) {
00092 temp.eval(dotstar, temp, x);
00093 dy.linComb(coef[i]*Scalar(i), temp);
00094 }
00095
00096 dy.eval(dotstar, dy,dx);
00097 }
00098
00100 virtual void applyAdjDeriv(const Vector<Scalar> & x,
00101 const Vector<Scalar> & dy,
00102 Vector<Scalar> & dx) const {
00103 applyDeriv(x,dy,dx);
00104 }
00105
00109 void applyInverseDeriv(const Vector<Scalar> & x,
00110 const Vector<Scalar> & dy,
00111 Vector<Scalar> & dx) const {
00112 Vector<Scalar> temp(x);
00113 ElementwiseMultiply<Scalar> dotstar;
00114 int size = coef.size();
00115 if( size > 1) {
00116 RVLAssignConst<Scalar> constantterm(coef[1]);
00117 dx.eval(constantterm);
00118 } else {
00119 dx.zero();
00120 }
00121 if( size > 2 )
00122 dx.linComb(coef[2]*Scalar(2), temp);
00123 for( int i = 3; i < size; i++) {
00124 temp.eval(dotstar, temp, x);
00125 dx.linComb(coef[i]*Scalar(i), temp);
00126 }
00127
00128 ElementwiseDivision<Scalar> dotslash;
00129 dx.eval(dotslash, dy, dx);
00130 }
00131
00133 void applyAdjInverseDeriv(const Vector<Scalar> & x,
00134 const Vector<Scalar> & dx,
00135 Vector<Scalar> & dy) const {
00136 applyInverseDeriv(x,dx,dy);
00137 }
00138
00139
00144 virtual Operator<Scalar> * clone() const { return new PolynomialOperator(*this); }
00145
00146 public:
00147
00148 PolynomialOperator(const std::valarray<Scalar> & _coef, Space<Scalar> & _spc)
00149 : spc(_spc), coef(_coef) {}
00150 PolynomialOperator(const PolynomialOperator<Scalar> & s): spc(s.spc), coef(s.coef) {}
00151 ~PolynomialOperator() {}
00152
00154 virtual const Space<Scalar> & getDomain() const { return spc; }
00155 virtual const Space<Scalar> & getRange() const { return spc; }
00156
00157 virtual void write(RVLException & e) const {
00158 e << "PolynomialOperator with coefficients";
00159 for(int i = 0; i < coef.size(); i++) {
00160 e << "c[" << i << "] = " << coef[i] << '\n';
00161 }
00162 }
00163
00164 virtual ostream & write(ostream & str) const {
00165 str << "PolynomialOperator with coefficients";
00166 for(int i = 0; i < coef.size(); i++) {
00167 str << "c[" << i << "] = " << coef[i] << '\n';
00168 }
00169 }
00170
00171
00172 };
00173
00174
00175 }
00176
00177 #endif