---

LUT.h

Go to the documentation of this file.

/* SimData: Data Infrastructure for Simulations
 * Copyright (C) 2002 Mark Rose <tm2@stm.lbl.gov>
 * 
 * This file is part of SimData.
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU 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 General Public License for more details.
 * 
 * You should have received a copy of the GNU 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.
 */


#include <SimData/PTS.h>

#ifdef __PTS_SIM__
        // old versions of msvc can't do partial template specialization
        #define __SIMDATA_NO_LUT__
        #define __SIMDATA_LUT_H__
#endif


#ifndef __SIMDATA_LUT_H__
#define __SIMDATA_LUT_H__


#include <vector>
#include <iostream>
#include <istream>
#include <cstdio>
#include <cmath>


#include <SimData/Archive.h>
#include <SimData/Exception.h>
#include <SimData/BaseType.h>


NAMESPACE_SIMDATA

SIMDATA_EXCEPTION(InterpolationInput);
SIMDATA_EXCEPTION(InterpolationError);
SIMDATA_EXCEPTION(InterpolationIndex);
SIMDATA_EXCEPTION(InterpolationUnpackMismatch);


// some simple debugging messages used during initial development
// LUTLOG() can be safely removed when no longer needed
#if 0
#  define __LUTLOG(msg) std::cout << msg << "\n";
#else
#  define __LUTLOG(msg)
#endif




// forward declaration for lookup table templates
template <int N, typename X>
class LUT;


template <int N, typename T>
class VEC {
        T m_Vec[N];
        T* m_Vp;
        int m_N;

        friend class VEC<N+1, T>;

        VEC(T* d): m_Vp(d), m_N(N) {}

public:
        VEC(std::vector<T> const &v): m_Vp(m_Vec), m_N(N) {
                if (v.size() != N) throw InterpolationIndex("VEC::VEC() dimension mismatch");
                for (int i = 0; i < N; i++) {
                        m_Vec[i] = v[i];
                }
        }

        VEC(T d): m_Vp(m_Vec), m_N(1) { m_Vec[0] = d; }

        VEC(): m_Vp(m_Vec), m_N(0) {}

        VEC<N,T> &set(T d) {
                m_N = 1;
                m_Vp = m_Vec;
                m_Vec[0] = d;
                return *this;
        }

        inline VEC<N,T> &operator()(T d) {
                if (m_N >= N) {
                        throw InterpolationIndex("VEC() too many coordinates specified");
                }
                m_Vec[m_N++] = d;
                return *this;
        }
        inline T &operator[](int n) const { 
                if (n < 0 || n >= m_N) {
                        throw InterpolationIndex("VEC[] index out of range");
                }
                return m_Vp[n]; 
        }

        inline VEC<N-1, T> rest() const {
                if (m_N != N) {
                        throw InterpolationIndex("VEC.rest() called for incomplete VEC");
                }
                return VEC<N-1, T>(m_Vp+1);
        }
};

template <typename T>
class VEC<1, T> {
friend class VEC<2, T>;
        T m_Vec;
        VEC(T *d): m_Vec(*d) {}
public:
        VEC(std::vector<T> const &v) {
                if (v.size() != 1) throw InterpolationIndex("VEC::VEC() dimension mismatch");
                m_Vec = v[0];
        }
        VEC(T d): m_Vec(d) {}
        inline void set(T d) { m_Vec = d; }
        inline T operator[](int n) const { 
                if (n != 0) {
                        throw InterpolationIndex("VEC[] index out of range");
                }
                return m_Vec; 
        }
};

template <typename T>
class VEC<0, T> {
        VEC(); // private ctor, not defined
};


template <int N, typename T>
class WRAP {
        LUT<N,T> const *m_Bind;
        VEC<N,T> m_Vec;
        int m_N;
        T m_Value;

friend class LUT<N,T>;

        WRAP(LUT<N,T> const *bind, T x): m_Bind(bind), m_Vec(x), m_N(1) {
                m_Value = static_cast<T>(0.0);
        }

public:
        inline WRAP<N,T> &operator[](T x);

        inline operator T() const { return m_Value; }
};


class SIMDATA_EXPORT Interpolation: public BaseType {
public:
        typedef enum { LINEAR, SPLINE } Modes;

        inline bool isInterpolated() const { return m_Interpolated; }

protected:
        Interpolation(): m_Interpolated(false) {}

        inline void checkInterpolated() const {
                if (!isInterpolated()) {
                        throw InterpolationError("LUT not interpolated");
                }
        }

        inline void checkNotInterpolated() const {
                if (isInterpolated()) {
                        throw InterpolationError("LUT already interpolated");
                }
        }

        void throwBreakpointOrder() const {
                throw InterpolationInput("Breakpoints must increase monotonically");
        }

        void throwInterpolationMode() const {
                throw InterpolationError("Unknown interpolation mode");
        }

        bool m_Interpolated;
};


template <typename X>
class SIMDATA_EXPORT InterpolationType: public Interpolation {
protected:
        X m_X0, m_X1, m_XS;
        int m_Limit;

        InterpolationType(): m_X0(0.0), m_X1(0.0), m_XS(0.0), m_Limit(1) {}

        inline void find(X x, int &i, X &f) const {
                x = (x - m_X0) * m_XS;
                i = static_cast<int>(floor(x)); 
                if (i >= m_Limit) {
                        i = m_Limit-1;
                        f = 1.0;
                } else
                if (i < 0) {
                        i = 0;
                        f = 0.0;
                } else {
                        f = x - i;
                }
        }

        void postInterpolation(X x0, X x1, int n);
};


template <int N, class X=float>
class Curvature;



template <int N, class X=float>
class SIMDATA_EXPORT LUT: public InterpolationType<X> {
        typedef LUT<N-1, X> Y;
        typedef std::vector< std::pair<X, Y> > DataVector;
        typedef std::vector<Y> TableVector;

        union {
                DataVector *m_Data;
                TableVector *m_Table;
        };

        int *m_Ref;

        inline std::pair<X, Y> &data(int n) {
                return (*m_Data)[n];
        }

        inline std::pair<X, Y> const &data(int n) const {
                return (*m_Data)[n];
        }

        inline int dataSize() const { return m_Data->size(); }

        inline Y &table(int n) {
                return (*m_Table)[n];
        }

        inline Y const &table(int n) const {
                return (*m_Table)[n];
        }

        inline int tableSize() const { return m_Table->size(); }

        inline void ref() {
                bool inc = false;
                if (isInterpolated()) {
                        inc = m_Table != NULL;
                } else {
                        inc = m_Data != NULL;
                }
                if (inc) {
                        assert(m_Ref);
                        *m_Ref += 1;
                }
        }

        inline void deref() {
                if (!m_Ref) return;
                *m_Ref -= 1;
                assert(*m_Ref >= 0);
                if (*m_Ref == 0) {
                        if (m_Interpolated) {
                                assert(m_Table);
                                __LUTLOG("~TABLE " << int(m_Table));
                                delete m_Table;
                                m_Table = 0;
                        } else {
                                assert(m_Data);
                                __LUTLOG("~DATA " << int(m_Data));
                                delete m_Data;
                                m_Data = 0;
                        }
                        delete m_Ref;
                        m_Ref = 0;
                }
        }

        inline void substitute(TableVector *_table) {
                assert(!isInterpolated() && _table);
                deref();
                m_Ref = new int(0);
                m_Table = _table;
                ref();
        }

        void createData(int n);
        void postInterpolation(TableVector *_table);
        void interpolateLinear(int n, TableVector *_table);
        void computeCurvature(Curvature<N,X> &c);
        void interpolateSpline(int n, TableVector *_table);

friend class LUT<N+1,X>;
private:
        // spline interpolation
        static void __splineInterpolate(X x, X h, LUT<N,X> const &y0, LUT<N,X> const &y1, 
                      Curvature<N,X> const &c0, 
                      Curvature<N,X> const &c1, 
                      LUT<N,X> &out);

        // 2 LUT parameter version does linear interpolation
        void __interpolate(X x, LUT<N,X> const &next, LUT<N,X> &out);
        
        X __getElement(int n);
        void __getDimension(std::vector<int> &dim);

public:
        typedef VEC<N, int> Dim;
        typedef VEC<N, X> Vec;
        typedef VEC<N, std::vector<X> > Breaks;

        LUT();

        LUT(LUT<N, X> const &copy);

        virtual ~LUT();

        LUT<N, X> const &operator=(LUT<N, X> const &copy);

        void interpolate(Dim const &dim, Interpolation::Modes mode);

        void interpolate(std::vector<int> const &dim, Interpolation::Modes mode);

        X getValue(Vec const &v) const;

        inline X getValue(std::vector<X> const &x) const {
                return getValue(Vec(x));
        }


        inline WRAP<N,X> operator[](X x) const {
                checkInterpolated();
                return WRAP<N,X>(this, x);
        }

        void load(std::vector<X> const &values, Breaks const &breaks, int *index = 0);
        void load(std::vector<X> const &values, std::vector< std::vector<X> > const &breaks);

        virtual void serialize(Archive&);

        virtual std::string asString() const;

        virtual std::string typeString() const;
};

template <typename X>
class SIMDATA_EXPORT LUT<1, X>: public InterpolationType<X> {
        typedef std::vector< std::pair<X, X> > DataVector;
        typedef std::vector<X> TableVector;

        union {
                DataVector *m_Data;
                TableVector *m_Table;
        };
        int *m_Ref;

        inline std::pair<X, X> &data(int n) {
                return (*m_Data)[n];
        }

        inline std::pair<X, X> const &data(int n) const {
                return (*m_Data)[n];
        }

        inline int dataSize() const { return m_Data->size(); }

        inline X &table(int n) {
                return (*m_Table)[n];
        }

        inline X const &table(int n) const {
                return (*m_Table)[n];
        }

        inline int tableSize() const { return m_Table->size(); }

        inline void ref() {
                if (m_Data) {
                        assert(m_Ref);
                        *m_Ref += 1;
                }
        }

        inline void deref() {
                if (!m_Ref) return;
                assert(m_Data);
                *m_Ref -= 1;
                assert(*m_Ref >= 0);
                if (*m_Ref == 0) {
                        if (m_Interpolated) {
                                assert(m_Table);
                                __LUTLOG("~TABLE " << int(m_Table));
                                delete m_Table;
                                m_Table = 0;
                        } else {
                                assert(m_Data);
                                __LUTLOG("~DATA " << int(m_Data));
                                delete m_Data;
                                m_Data = 0;
                        }
                        delete m_Ref;
                        m_Ref = 0;
                }
        }

        inline void substitute(TableVector *_table) {
                assert(!isInterpolated() && _table);
                deref();
                m_Ref = new int(0);
                m_Table = _table;
                ref();
        }

        void createData(int n);
        void postInterpolation(TableVector *_table);
        void interpolateLinear(int n, TableVector *_table);
        void interpolateSpline(int n, TableVector *_table);

friend class LUT<2,X>;
private:

        inline X __getElement(int n) {
                return data(n).second;
        }

        inline void __getDimension(std::vector<int> &dim) {
                dim.push_back(dataSize());
        }

        // linear interpolation between this and next
        void __interpolate(X x, LUT<1,X> const &next, LUT<1,X> &out);

        // spline interpolation 
        static void __splineInterpolate(
                        X x, X h, 
                        LUT<1,X> const &y0, 
                        LUT<1,X> const &y1, 
                        Curvature<1,X> const &c0, 
                        Curvature<1,X> const &c1, 
                        LUT<1,X> &out);

public:
        typedef VEC<1, int> Dim;
        typedef VEC<1, X> Vec;
        typedef VEC<1, std::vector<X> > Breaks;

        inline X __get(int n) const { return data(n).second; }

        LUT();
        virtual ~LUT();
        LUT(LUT<1, X> const &copy);
        LUT<1, X> const &operator=(LUT<1, X> const &copy);

        void interpolate(Dim const &dim, Interpolation::Modes mode);
        void interpolate(std::vector<int> const &dim, Interpolation::Modes mode);

        X getValue(Vec const &v) const;

        inline X getValue(std::vector<X> const &x) const {
                return getValue(Vec(x));
        }

        inline X operator[](Vec const &v) const {
                return getValue(v);
        }

        void load(std::vector<X> const &values, Breaks const &breaks, int *index = 0);
        void load(std::vector<X> const &values, std::vector< std::vector<X> > const &breaks);

        virtual void serialize(Archive&);

        virtual std::string asString() const;

        virtual std::string typeString() const;
};

template <typename X>
class SIMDATA_EXPORT LUT<0, X>: public Interpolation {
        LUT();
};


template <int N, typename X>
inline WRAP<N,X> &WRAP<N,X>::operator[](X x) {
        m_Vec(x);
        if (++m_N == N && m_Bind) {
                m_Value = m_Bind->getValue(m_Vec);
        }
        return *this;
}


typedef LUT<1,float> Table1;

typedef LUT<2,float> Table2;

typedef LUT<3,float> Table3;


NAMESPACE_SIMDATA_END // simdata


#endif // __SIMDATA_LUT_H__

---

SimData version pre-0.4.0. For more information on SimData, visit the SimData Homepage. Generated on Tue Oct 14 12:06:38 2003, using Doxygen 1.2.18.