Paraboloid.h

Go to the documentation of this file.
00001 //===========================================================================
00040 //===========================================================================
00041 
00042 //=============================================================
00043 //
00044 // paraboloid test function
00045 //
00046 
00047 #ifndef PARABOLOID_H
00048 #define PARABOLOID_H
00049 
00050 #include <SharkDefs.h>
00051 #include <ReClaM/Model.h>
00052 #include <ReClaM/ErrorFunction.h>
00053 #include <Rng/GlobalRng.h>
00054 #include <vector>
00055 #include <Array/ArrayOp.h>
00056 #include <Array/ArrayIo.h>
00057 
00058 
00066 class Paraboloid : public Model, public ErrorFunction
00067 {
00068 public:
00069     Paraboloid(unsigned _d, double _c = 10, bool basis = true)
00070     {
00071         unsigned i;
00072         d = _d;
00073         cond = _c;
00074 
00075         parameter.resize(d);
00076         parameter = 0.0;
00077         B   .resize(d, d);
00078         if (basis) generateBasis();
00079         else
00080         {
00081             B  = 0;
00082             for (i = 0; i < d; i++) B(i, i) = 1;
00083         }
00084     }
00085 
00086     void model(const Array<double> &input, Array<double> &output)
00087     {
00088     }
00089 
00090     void modelDerivative(const Array<double> &input, Array<double>& derivative)
00091     {
00092     }
00093 
00094     double error(const std::vector<double> &v)
00095     {
00096         unsigned i;
00097         if (v.size() != d)
00098         {
00099             throw SHARKEXCEPTION("dimension mismatch");
00100         }
00101         double sum = 0.;
00102         for (i = 0; i < d; i++) parameter(i) = v[i];
00103         for (i = 0; i < d; i++)
00104         {
00105             sum += Shark::sqr(pow(cond, (double(i) / double(d - 1))) * scalarProduct(parameter, B.col(i))) ;
00106         }
00107         return sum;
00108     }
00109 
00110     double error(Model& model, const Array<double> &input, const Array<double> &target)
00111     {
00112         double sum = 0.;
00113         unsigned i;
00114         for (i = 0; i < d; i++)
00115             sum += Shark::sqr(pow(cond, double(i) / double(d - 1)) * scalarProduct(parameter, B.col(i)));
00116         return sum;
00117     }
00118 
00119     double errorDerivative(Model& model, const Array<double> &input, const Array<double> &target, Array<double>& derivative)
00120     {
00121         double sum = 0.;
00122         unsigned k, i;
00123 
00124         derivative.resize(d, false);
00125         derivative = 0.0;
00126         for (k = 0; k < d; k++)
00127         {
00128             sum += Shark::sqr(pow(1.0, (double) k / (double)(d - 1)) * scalarProduct(parameter, B.col(k)));
00129             for (i = 0; i < d; i++)
00130                 derivative(k) += 2.0 * pow(cond, 2.0 * double(i) / double(d - 1)) * scalarProduct(parameter, B.col(i)) * B(k, i);
00131         }
00132         return sum;
00133     }
00134 
00135     void init(unsigned s, double a, double b, bool basis = true)
00136     {
00137         unsigned i;
00138         Rng::seed(s);
00139         if (basis) generateBasis();
00140         for (i = 0; i < d; i++) parameter(i) = Rng::uni(a, b);
00141     }
00142 
00143     Array<double> &getBasis()
00144     {
00145         return B;
00146     }
00147 
00148 protected:
00149     void generateBasis()
00150     {
00151         unsigned i, j, c;
00152         Array<double> H;
00153         B.resize(d, d);
00154         H.resize(d, d);
00155         for (i = 0; i < d; i++)
00156         {
00157             for (c = 0; c < d; c++)
00158             {
00159                 H(i, c) = Rng::gauss(0, 1);
00160                 //H(i, c) = (c==i) ? 1 : 0;
00161             }
00162         }
00163         B = H;
00164         for (i = 0; i < d; i++)
00165         {
00166             for (j = 0; j < i; j++)
00167                 for (c = 0; c < d; c++)
00168                     B(i, c) -= scalarProduct(H[i], H[j]) * H(j, c) / scalarProduct(H[j], H[j]);
00169             H = B;
00170         }
00171         for (i = 0; i < d; i++)
00172         {
00173             double normB = sqrt(scalarProduct(B[i], B[i]));
00174             for (j = 0; j < d; j++)
00175                 B(i, j) = B(i, j) / normB;
00176         }
00177     }
00178 
00179     Array<double > B;
00180     unsigned       d;
00181     double         cond;
00182 };
00183 
00184 #endif
00185 
00186 
00187 
00188