STK_LinAlgebra1D.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2007  Serge Iovleff

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this program; if not, write to the
    Free Software Foundation, Inc.,
    59 Temple Place,
    Suite 330,
    Boston, MA 02111-1307
    USA

    Contact : Serge.Iovleff@stkpp.org
*/

/*
 * Project:  stkpp::Algebra
 * Purpose:  Define 1D Linear Algebra methods with Real Containers
 * Author:   Serge Iovleff, serge.iovleff@stkpp.org
 *
 **/

#ifndef STK_LINALGEBRA1D_H
#define STK_LINALGEBRA1D_H

#include "../../STKernel/include/STK_Real.h"
#include "../../STKernel/include/STK_Misc.h"

#include "../../Sdk/include/STK_ITContainer1D.h"

namespace STK
{

template<class Container1D>
Real sum( ITContainer1D<Real, Container1D> const& x)
{
  // get dimensions
  const Integer first = x.first(), last = x.last();
  // compute the sum
  Real sum = 0.0;
  for (Integer i=first; i<=last; i++)
    sum += x[i];
  return (sum);
}

template<class Container1D1, class Container1D2>
Real weightedSum( ITContainer1D<Real, Container1D1> const& x
                , ITContainer1D<Real, Container1D2> const& w
                )
{
#ifdef STK_DEBUG
    if (!x.range().isIn(w.range()))
      throw runtime_error("In weightedSum(x, w) "
                               "x.range() not include in w.range()");
#endif
  // dimensions
  const Integer first = x.first(), last = x.last();
  // compute the weighted sum
  Real sum = 0.0;
  for (Integer i=first; i<=last; i++)
    sum += w[i] * x[i];
  return (sum);
}

template<class Container1D>
Real normInf( ITContainer1D<Real, Container1D> const& x)
{
  // get dimensions
  const Integer first = x.first(), last = x.last();
  // compute the maxmal value
  Real scale = 0.0;
  for (Integer i=first; i<=last; i++)
    scale = max(scale, abs(x[i]));
  return (scale);
}

template<class Container1D1, class Container1D2>
Real weightedNormInf( ITContainer1D<Real, Container1D1> const& x
                    , ITContainer1D<Real, Container1D2> const& w
                    )
{
#ifdef STK_DEBUG
    if (!x.range().isIn(w.range()))
      throw runtime_error("In weightedNormInf(x, w) "
                               "x.range() not include in w.range()");
#endif
  // get dimensions
  const Integer first = x.first(), last = x.last();
  // compute weighted norm inf
  Real scale = 0.0;
  for (Integer i=first; i<=last; i++)
    scale = max(scale, abs(w[i]*x[i]));
  return (scale);
}

template<class Container1D>
Real normTwo( ITContainer1D<Real, Container1D> const& x)
{
  // compute the maximal value of x
  Real scale =normInf(x), norm =0.0;
  if (scale)
  {
    // get dimensions
    const Integer first = x.first(), last = x.last();
    // sum squared normalized values
    for (Integer i = first; i<=last; i++)
    {
      const Real aux = x[i]/scale;
      norm += aux * aux;
    }
  }
  // rescale sum
  return (Real(sqrt(double(norm)))*scale);
}

template<class Container1D1, class Container1D2>
Real weightedNormTwo( ITContainer1D<Real, Container1D1> const& x
                    , ITContainer1D<Real, Container1D2> const& w
                    )
{
#ifdef STK_DEBUG
    if (!x.range().isIn(w.range()))
      throw runtime_error("In weightedNormTwo(x, w) "
                               "x.range() not include in w.range()");
#endif
  // compute the maximal value of x
  Real scale = weightedNormInf(x, w), norm2 =0.0;
  if (scale)
  {
    // get dimensions
    const Integer first = x.first(), last = x.last();
    // compute norm2
    for (Integer i = first; i<=last; i++)
    {
      const Real aux = (w[i]*x[i])/scale;
      norm2 += aux * aux;
    }
  }
  // rescale sum
  return (Real(sqrt(double(norm2)))*scale);
}

template<class Container1D>
Real normTwo2( ITContainer1D<Real, Container1D> const& x)
{
  Real scale =normInf(x), norm =0.0;
  if (scale)
  {
    // get dimensions
    const Integer first = x.first(), last = x.last();
    // sum squared normalized values
    for (Integer i = first; i<=last; i++)
    {
      Real aux = x[i]/scale;
      norm += aux * aux;
    }
  }
  // scale result
  return (norm*scale*scale);
}

template<class Container1D1, class Container1D2>
Real weightedNormTwo2( ITContainer1D<Real, Container1D1> const& x
                     , ITContainer1D<Real, Container1D2> const& w
                     )
{
#ifdef STK_DEBUG
  if (!x.range().isIn(w.range()))
    throw runtime_error("In weightedNormTwo2(x, w) "
                             "x.range() not include in w.range()");
#endif
  Real scale =weightedNormInf(x, w), norm =0.0;
  if (scale)
  {
    // get dimensions
    const Integer first = x.first(), last = x.last();
    // compute norm2
    for (Integer i = first; i <= last; i++)
    {
      const Real aux = x[i]/scale;
      norm += w[i] * aux * aux;
    }
  }
  // scale result
  return (norm*scale*scale);
}

template<class Container1D1, class Container1D2>
Real dot( ITContainer1D<Real, Container1D1> const& x
        , ITContainer1D<Real, Container1D2> const& y
        )
{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
  // compute the sum product
  Real sum=0.0;
  Integer i;
  for (i = first; i<last; i+=2)
    sum += x[i] * y[i] + x[i+1] * y[i+1];
  // check if the number of element is odd
  if (i==last) sum+=x[last]*y[last];
  return (sum);
}

template<class Container1D1, class Container1D2, class Container1D3>
Real weightedDot( ITContainer1D<Real, Container1D1> const& x
                , ITContainer1D<Real, Container1D2> const& y
                , ITContainer1D<Real, Container1D3> const& w
                )
{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
#ifdef STK_DEBUG
  if (!Range(first,last).isIn(w.range()))
    throw runtime_error("In weightedDot(x, w) "
                             "Range(first,last) not include in w.range()");
#endif
  // compute the sum product
  Real sum=0.0;
  Integer i;
  for (i = first; i<last; i+=2)
    sum += w[i]*x[i] * y[i] + w[i+1]*x[i+1] * y[i+1];
  // check if last is odd
  if (i == last) sum += w[last]*x[last]*y[last];
  return (sum);
}

template<class Container1D1, class Container1D2>
Real dist( ITContainer1D<Real, Container1D1> const& x
         , ITContainer1D<Real, Container1D2> const& y
         )
{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
  // compute the maximal difference
  Real scale = 0.;
  for (Integer i = first; i<=last; i++)
    scale = max(scale, abs(x[i] - y[i]));
  // Compute the norm
  Real norm2 = 0.;
  if (scale)
  { // comp the norm^2
    for (Integer i = first; i<=last; i++)
    {
      const Real aux = (x[i]-y[i])/scale;
      norm2 += aux * aux;
    }
  }
  // rescale sum
  return (Real(sqrt(double(norm2)))*scale);
}

template<class Container1D1, class Container1D2, class Container1D3>
Real weightedDist( ITContainer1D<Real, Container1D1> const& x
                 , ITContainer1D<Real, Container1D2> const& y
                 , ITContainer1D<Real, Container1D3> const& w
                 )
{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
#ifdef STK_DEBUG
  if (!Range(first,last).isIn(w.range()))
    throw runtime_error("In weightedDist(x, w) "
                             "Range(first,last) not include in w.range()");
#endif
  // compute the maximal difference
  Real scale = 0., norm2= 0.;
  for (Integer i = first; i<=last; i++)
    scale = max(scale, abs(w[i]*(x[i] - y[i])));
  // Compute the norm
  if (scale)
  { // comp the norm^2
    for (Integer i = first; i<=last; i++)
    {
      const Real aux = (x[i]-y[i])/scale;
      norm2 += w[i]*aux * aux;
    }
  }
  // rescale sum
  return (Real(sqrt(double(norm2)))*scale);
}

template<class Container1D1, class Container1D2>
void add(ITContainer1D<Real, Container1D1> const& x, ITContainer1D<Real, Container1D2>& y)
{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
  for (Integer i = first; i<=last; i++)
  {
    y[i] += x[i];
  }
}

template<class Container1D1, class Container1D2, class Container1D3>
void add( ITContainer1D<Real, Container1D1> const& x
        , ITContainer1D<Real, Container1D2> const& y
        , ITContainer1D<Real, Container1D3>& z
        )
{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
  z.resize(Range(first, last));
  for (Integer i = first; i<=last; i++)
  {
    z[i] = x[i] + y[i];
  }
}

template<class Container1D1, class Container1D2>
void diff(ITContainer1D<Real, Container1D1> const& x, ITContainer1D<Real, Container1D2>& y)
{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
  for (Integer i = first; i<=last; i++)
  {
    y[i] -= x[i];
  }
}

template<class Container1D1, class Container1D2, class Container1D3>
void diff( ITContainer1D<Real, Container1D1> const& x
         , ITContainer1D<Real, Container1D2> const& y
         , ITContainer1D<Real, Container1D3>& z
         )
{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
  z.resize(Range(first, last));
  for (Integer i = first; i<=last; i++)
  {
    z[i] = x[i] - y[i];
  }
}


} // Namespace STK

#endif // STK_LINALGEBRA1D_H