PLearn 0.1
CubicSpline.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // CubicSpline.cc
00004 //
00005 // Copyright (C) 2007 Christian Dorion
00006 //
00007 // Redistribution and use in source and binary forms, with or without
00008 // modification, are permitted provided that the following conditions are met:
00009 //
00010 //  1. Redistributions of source code must retain the above copyright
00011 //     notice, this list of conditions and the following disclaimer.
00012 //
00013 //  2. Redistributions in binary form must reproduce the above copyright
00014 //     notice, this list of conditions and the following disclaimer in the
00015 //     documentation and/or other materials provided with the distribution.
00016 //
00017 //  3. The name of the authors may not be used to endorse or promote
00018 //     products derived from this software without specific prior written
00019 //     permission.
00020 //
00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031 //
00032 // This file is part of the PLearn library. For more information on the PLearn
00033 // library, go to the PLearn Web site at www.plearn.org
00034 
00035 // Authors: Christian Dorion
00036 
00040 #include "CubicSpline.h"
00041 //#include "algo.h"
00042 
00043 namespace PLearn {
00044 using namespace std;
00045 
00046 PLEARN_IMPLEMENT_OBJECT(
00047     CubicSpline,
00048     "Unidimensional cubic spline learner.",
00049     "This learner fits a unidimensional cubic spline to a given set of\n"
00050     "points. That is, inputsize() must be one and the inputs are considered to be\n"
00051     "the x values, while targetsize() must also be one and the targets are\n"
00052     "considered to be the y values. The spline is fitted to the (x,y)-pairs so\n"
00053     "formed. X values don't need to be ordered; the ordering is ensured within\n"
00054     "the train method.\n");
00055 
00056 CubicSpline::CubicSpline()
00057     : m_low_slope(MISSING_VALUE),
00058       m_high_slope(MISSING_VALUE)
00059 {
00060     // Nothing to do here
00061 } 
00062 
00063 void CubicSpline::declareOptions(OptionList& ol)
00064 {
00065     declareOption(ol, "low_slope", &CubicSpline::m_low_slope,
00066                   OptionBase::buildoption,
00067                   "The slope to enforce at the leftmost node -- Default: NaN [None]");
00068     
00069     declareOption(ol, "high_slope", &CubicSpline::m_high_slope,
00070                   OptionBase::buildoption,
00071                   "The slope to enforce at the rightmost node -- Default: NaN [None]");
00072     
00073     // Learnt options
00074     declareOption(ol, "inputs", &CubicSpline::m_inputs,
00075                   OptionBase::learntoption,
00076                   "A buffer containing the last training set inputs");
00077 
00078     declareOption(ol, "targets", &CubicSpline::m_targets,
00079                   OptionBase::learntoption,
00080                   "A buffer containing the last training set targets");
00081 
00082     declareOption(ol, "coefficients", &CubicSpline::m_coefficients,
00083                   OptionBase::learntoption,
00084                   "The learnt coefficients");
00085     
00086     // Now call the parent class' declareOptions
00087     inherited::declareOptions(ol);
00088 }
00089 
00090 void CubicSpline::build_()
00091 {
00092     // Nothing to do here    
00093 }
00094 
00095 // ### Nothing to add here, simply calls build_
00096 void CubicSpline::build()
00097 {
00098     inherited::build();
00099     build_();
00100 }
00101 
00102 
00103 void CubicSpline::makeDeepCopyFromShallowCopy(CopiesMap& copies)
00104 {
00105     inherited::makeDeepCopyFromShallowCopy(copies);
00106 
00107     deepCopyField(m_inputs, copies);
00108     deepCopyField(m_targets, copies);
00109     deepCopyField(m_coefficients, copies);
00110 }
00111 
00112 
00113 int CubicSpline::outputsize() const
00114 {
00115     return inputsize();
00116 }
00117 
00118 void CubicSpline::forget()
00119 {
00120     m_coefficients->resize(0);
00121     inherited::forget();
00122 }
00123 
00124 void CubicSpline::train()
00125 {
00126     // This learner fits unidimensional inputs/targets
00127     PLASSERT( inputsize() == 1 );
00128     PLASSERT( targetsize() == 1 );
00129     
00130     // Train set is a member of PLearner; set through setTrainingSet()
00131     int n = train_set->length();
00132     m_inputs = train_set.getColumn(0);
00133     m_targets = train_set.getColumn(1);
00134     PLASSERT( n >= 2 && train_set->width() == 2 );
00135 
00136     // Sort the inputs and targets along the inputs values
00137     TVec<int> indices = m_inputs.sortingPermutation();
00138     m_inputs          = m_inputs(indices);
00139     m_targets         = m_targets(indices);
00140     
00141     Vec u(n-1, 0.0);
00142     m_coefficients.resize(n);
00143     m_coefficients.fill(0.0);
00144 
00145     // Low slope hack
00146     if ( !is_missing(m_low_slope) ) {
00147         u[0] = (3.0/(m_inputs[1]-m_inputs[0])) *
00148             ( (m_targets[1]-m_targets[0])/
00149               (m_inputs[1] - m_inputs[0]) - m_low_slope );
00150         m_coefficients[0] = -0.5;
00151     }
00152 
00153     // Forward pass on coefficients
00154     for (int i=1; i < (n-1); i++) {
00155         real sig = (m_inputs[i]-m_inputs[i-1]) / (m_inputs[i+1]-m_inputs[i-1]);        
00156         real p   = sig*m_coefficients[i-1]+2.0;
00157         
00158         m_coefficients[i] = (sig-1.0)/p;
00159         
00160         u[i] = (m_targets[i+1]-m_targets[i]) / (m_inputs[i+1]-m_inputs[i])
00161             - (m_targets[i]-m_targets[i-1])  / (m_inputs[i]-m_inputs[i-1]);
00162         
00163         u[i] = (6.0*u[i]/(m_inputs[i+1]-m_inputs[i-1]) - sig*u[i-1]) / p;
00164     }
00165 
00166     // High slope hack
00167     real un=0, qn=0;
00168     if ( !is_missing(m_high_slope) ) {
00169         qn = 0.5;
00170         un = (3.0/(m_inputs[n-1]-m_inputs[n-2])) *
00171             (m_high_slope  -  (m_targets[n-1]-m_targets[n-2]) / (m_inputs[n-1]-m_inputs[n-2]));
00172     }
00173 
00174     // Compute the last coefficient
00175     m_coefficients[n-1] = (un-qn*u[n-2])/(qn*m_coefficients[n-1]+1.0);
00176 
00177     // Backsubstitution step
00178     for (int k=n-2; k >= 0; k--)
00179         m_coefficients[k] = m_coefficients[k]*m_coefficients[k+1]+u[k];
00180 }
00181 
00182 
00183 void CubicSpline::computeOutput(const Vec& input, Vec& output) const
00184 {
00185     PLASSERT( input.length() == 1 );
00186     output.resize(1);
00187  
00188     int n = m_inputs.length();
00189     real x = input[0];
00190     int pos = min( max(1, m_inputs.findSorted(x)),  n-1 );
00191 
00192     real h = m_inputs[pos] - m_inputs[pos-1];
00193     PLASSERT( h > 0.0 );
00194 
00195     real a = (m_inputs[pos] - x) / h;
00196     real b = (x - m_inputs[pos-1]) / h;
00197     output[0] = a*m_targets[pos-1] + b*m_targets[pos]
00198         + ((a*a*a - a)*m_coefficients[pos-1]
00199            + (b*b*b - b)*m_coefficients[pos]) * h*h/6.0;
00200 }
00201 
00202 void CubicSpline::computeCostsFromOutputs(const Vec& input, const Vec& output,
00203                                            const Vec& target, Vec& costs) const
00204 {
00205     // No costs...
00206 }
00207 
00208 TVec<string> CubicSpline::getTestCostNames() const
00209 {
00210     // None for now
00211     return TVec<string>();
00212 }
00213 
00214 TVec<string> CubicSpline::getTrainCostNames() const
00215 {
00216     // None for now
00217     return TVec<string>();
00218 }
00219 
00220 
00221 } // end of namespace PLearn
00222 
00223 
00224 /*
00225   Local Variables:
00226   mode:c++
00227   c-basic-offset:4
00228   c-file-style:"stroustrup"
00229   c-file-offsets:((innamespace . 0)(inline-open . 0))
00230   indent-tabs-mode:nil
00231   fill-column:79
00232   End:
00233 */
00234 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines