polyop.hh

Go to the documentation of this file.
00001 /*************************************************************************
00002 
00003 Copyright Rice University, 2004.
00004 All rights reserved.
00005 
00006 Permission is hereby granted, free of charge, to any person obtaining a
00007 copy of this software and associated documentation files (the "Software"),
00008 to deal in the Software without restriction, including without limitation
00009 the rights to use, copy, modify, merge, publish, distribute, and/or sell
00010 copies of the Software, and to permit persons to whom the Software is
00011 furnished to do so, provided that the above copyright notice(s) and this
00012 permission notice appear in all copies of the Software and that both the
00013 above copyright notice(s) and this permission notice appear in supporting
00014 documentation.
00015 
00016 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
00017 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
00018 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY
00019 RIGHTS. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS
00020 NOTICE BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL
00021 DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
00022 PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
00023 ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
00024 THIS SOFTWARE.
00025 
00026 Except as contained in this notice, the name of a copyright holder shall
00027 not be used in advertising or otherwise to promote the sale, use or other
00028 dealings in this Software without prior written authorization of the
00029 copyright holder.
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       // dy has DF(x), so we need to do the multiply now.
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       // dy has DF(x), so we need to do the division now.
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

Generated on 5 Jan 2017 for LocalRVL by  doxygen 1.4.7