gentensor.h

Go to the documentation of this file.

/*
  This file is part of MADNESS.

  Copyright (C) 2007,2010 Oak Ridge National Laboratory

  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

  For more information please contact:

  Robert J. Harrison
  Oak Ridge National Laboratory
  One Bethel Valley Road
  P.O. Box 2008, MS-6367

  email: harrisonrj@ornl.gov
  tel:   865-241-3937
  fax:   865-572-0680

  $Id$
*/

#ifndef GENTENSOR_H_
#define GENTENSOR_H_


#include <madness/tensor/tensor.h>
#include <madness/tensor/srconf.h>
#include <madness/tensor/tensortrain.h>
#include <stdexcept>

namespace madness {

        // a forward declaration
        template <class T> class SliceGenTensor;
        template <class T> class GenTensor;


        struct TensorArgs {
                double thresh;
                TensorType tt;
        TensorArgs() : thresh(-1.0), tt(TT_NONE) {}
                TensorArgs(const double& thresh1, const TensorType& tt1)
                        : thresh(thresh1)
                        , tt(tt1) {
                }
                static std::string what_am_i(const TensorType& tt) {
                        if (tt==TT_2D) return "TT_2D";
                        if (tt==TT_FULL) return "TT_FULL";
                        return "unknown tensor type";
                }
                template <typename Archive>
                void serialize(const Archive& ar) {
                    int i=int(tt);
                    ar & thresh & i;
                    tt=TensorType(i);
                }
        };


    template<typename T>
    bool has_zero_rank(const GenTensor<T>& g) {
        return ((g.rank()==0) or (g.tensor_type()==TT_NONE));
    }


#if !HAVE_GENTENSOR

        template <typename T>
        class GenTensor : public Tensor<T> {

        public:

                GenTensor<T>() : Tensor<T>() {}

                GenTensor<T>(const Tensor<T>& t1) : Tensor<T>(t1) {}
                GenTensor<T>(const Tensor<T>& t1, const TensorArgs& targs) : Tensor<T>(t1) {}
                GenTensor<T>(const Tensor<T>& t1, double eps, const TensorType tt) : Tensor<T>(t1) {}
                GenTensor<T>(const TensorType tt): Tensor<T>() {}
                GenTensor<T>(std::vector<long> v, const TensorType& tt) : Tensor<T>(v) {}
                GenTensor<T>(std::vector<long> v, const TensorArgs& targs) : Tensor<T>(v) {}
                GenTensor<T>(const SRConf<T>& sr1) : Tensor<T>() {MADNESS_EXCEPTION("no ctor with SRConf: use HAVE_GENTENSOR",1);}

        template <class Q> operator GenTensor<Q>() const { // type conv => deep copy
            Tensor<Q> result = Tensor<Q>(this->_ndim,this->_dim,false);
            BINARY_OPTIMIZED_ITERATOR(Q, result, const T, (*this), *_p0 = (Q)(*_p1));
            return result;
        }


                GenTensor<T> reconstruct_tensor() const {return *this;}
                GenTensor<T> full_tensor() const {return *this;}
                GenTensor<T>& full_tensor() {return *this;}
        GenTensor<T> full_tensor_copy() const {return *this;}
        GenTensor<T> full_tensor_copy() {return *this;}

        bool has_data() const {return this->size()>0;};
        bool has_no_data() const {return not has_data();};
                long rank() const {return -1;}
                double svd_normf() const {return this->normf();}
                size_t real_size() const {return this->size();}

        void reduce_rank(const double& eps) {return;};
        void normalize() {return;}

        std::string what_am_i() const {return "GenTensor, aliased to Tensor";};
                TensorType tensor_type() const {return TT_FULL;}

                void add_SVD(const GenTensor<T>& rhs, const double& eps) {*this+=rhs;}

                SRConf<T> config() const {MADNESS_EXCEPTION("no SRConf in complex GenTensor",1);}
        SRConf<T> get_configs(const int& start, const int& end) const {MADNESS_EXCEPTION("no SRConf in complex GenTensor",1);}

        template<typename Q>
        GenTensor<T> general_transform(const Tensor<Q> c[]) const {

            return madness::general_transform(*static_cast<const Tensor<T>*>(this),c);
        }


                static double fac_reduce() {return -1.0;};

        };

    template <class T>
        GenTensor<T> reduce(std::list<GenTensor<T> >& addends, double eps, bool are_optimal=false) {
        typedef typename std::list<GenTensor<T> >::iterator iterT;
        GenTensor<T> result=copy(addends.front());
        for (iterT it=++addends.begin(); it!=addends.end(); ++it) {
                result+=*it;
        }
            return result;
    }


    template <class T>
    GenTensor<T> outer(const GenTensor<T>& left, const GenTensor<T>& right) {
        return madness::outer(static_cast<Tensor<T> >(left),static_cast<Tensor<T> >(right));
    }


     template <class T>
     GenTensor<T> outer_low_rank(const Tensor<T>& left, const Tensor<T>& right) {
         return madness::outer(left,right);
     }


        namespace archive {
        template <class Archive, typename T>
        struct ArchiveStoreImpl< Archive, GenTensor<T> > {
                static void store(const Archive& s, const GenTensor<T>& t) {
                        if (t.iscontiguous()) {
                                s & t.size() & t.id();
                                if (t.size()) s & t.ndim() & wrap(t.dims(),TENSOR_MAXDIM) & wrap(t.ptr(),t.size());
                        }
                        else {
                                s & copy(t);
                        }
                };
        };


        template <class Archive, typename T>
        struct ArchiveLoadImpl< Archive, GenTensor<T> > {
                static void load(const Archive& s, GenTensor<T>& t) {
                        long sz = 0l, id =0l;
                        s & sz & id;
                        if (id != t.id()) throw "type mismatch deserializing a tensor";
                        if (sz) {
                                long _ndim=0l, _dim[TENSOR_MAXDIM];
                                s & _ndim & wrap(_dim,TENSOR_MAXDIM);
                                t = Tensor<T>(_ndim, _dim, false);
                                if (sz != t.size()) throw "size mismatch deserializing a tensor";
                                s & wrap(t.ptr(), t.size());
                        }
                        else {
                                t = Tensor<T>();
                        }
                };
        };

        }

#else



        template<typename T>
        class GenTensor {

        private:

                friend class SliceGenTensor<T>;

                typedef SRConf<T> configT;
                typedef Tensor<T> tensorT;
                typedef GenTensor<T> gentensorT;
                typedef std::shared_ptr<configT> sr_ptr;

                sr_ptr _ptr;

                static double machinePrecision() {return 1.e-14;}

                static double facReduce() {return 1.e-3;}

        public:

                GenTensor() : _ptr() {
                }

//              GenTensor(const GenTensor<T>& rhs) : _ptr(rhs._ptr) { // DON'T DO THIS: USE_COUNT BLOWS UP
                GenTensor(const GenTensor<T>& rhs) : _ptr() {
                        if (rhs.has_data()) _ptr=rhs._ptr;
                };

                GenTensor(const std::vector<long>& dim, const TensorType tt) : _ptr() {

                // check consistency
                const long ndim=dim.size();
                const long maxk=dim[0];
                for (long idim=0; idim<ndim; idim++) {
                        MADNESS_ASSERT(maxk==dim[0]);
                }

                        _ptr=sr_ptr(new configT(dim.size(),dim[0],tt));

                }

                GenTensor(const std::vector<long>& dim, const TensorArgs& targs) : _ptr() {

                        // check consistency
                const long ndim=dim.size();
                const long maxk=dim[0];
                for (long idim=0; idim<ndim; idim++) {
                        MADNESS_ASSERT(maxk==dim[0]);
                }

                        _ptr=sr_ptr(new configT(dim.size(),dim[0],targs.tt));

                }

                GenTensor(const TensorType& tt, const unsigned int& k, const unsigned int& dim) {
                        _ptr=sr_ptr(new configT(dim,k,tt));
                }

                GenTensor(const Tensor<T>& rhs, const TensorArgs& targs) {

                        // fast return if possible
                        if (not rhs.has_data()) {
                                _ptr.reset();
                                return;
                        }

                        MADNESS_ASSERT(rhs.ndim()>0);
                        for (long idim=0; idim<rhs.ndim(); idim++) {
                            MADNESS_ASSERT(rhs.dim(0)==rhs.dim(idim));
                        }

                        _ptr=sr_ptr(new configT(rhs.ndim(),rhs.dim(0),targs.tt));

                        // direct reduction on the polynomial values on the Tensor
                        TensorType ttype=tensor_type();
                        if (ttype==TT_2D) {

#if 1
                                Tensor<T> U,VT;
                                Tensor< typename Tensor<T>::scalar_type > s;

                                // add some extra dimensions if this is supposedly NS form:
                                // separate sum and wavelet coeffs
                                if (rhs.dim(0)%2==0) {
                                        std::vector<long> dims(rhs.ndim()*2);
                                        for (int i=0; i<rhs.ndim(); ++i) {
                                                int k=rhs.dim(i);
                                                dims[2*i]=k/2;
                                                dims[2*i+1]=2;
                                        }
                                        TensorTrain<T> tt(rhs,targs.thresh*facReduce(),dims);
                                        // fuse sum and wavelet coeffs back together
                                        for (int i=0; i<rhs.ndim(); ++i) tt.fusedim(i);
                                        // fuse dimensions into particles 1 and 2
                                        tt.two_mode_representation(U,VT,s);

                                } else {
                                        TensorTrain<T> tt(rhs,targs.thresh*facReduce());
                                        tt.two_mode_representation(U,VT,s);
                                }

                                const long r=VT.dim(0);
                                const long nd=VT.ndim();
                                if (r==0) {
                    _ptr=sr_ptr(new configT(dim(),get_k(),ttype));
                                } else {
                                        MADNESS_ASSERT(U.dim(nd-1)==r);
                                        Tensor<T> UU=U.reshape(U.size()/r,r);
                                        _ptr=sr_ptr(new configT(s, copy(transpose(UU)), VT.reshape(r,VT.size()/r),
                                                          dim(), get_k()));
                                        this->normalize();
                                }
#else
                                // adapt form of values
                                MADNESS_ASSERT(rhs.iscontiguous());
                                std::vector<long> d(_ptr->dim_eff(),_ptr->kVec());
                                Tensor<T> values_eff=rhs.reshape(d);

                                this->computeSVD(targs.thresh,values_eff);
#endif
                        } else if (ttype==TT_FULL){
                                _ptr.reset(new configT(copy(rhs)));
                        } else {
                                MADNESS_EXCEPTION("unknown tensor type in GenTensor(tensorT,targs)",0);
                        }
                }

                GenTensor(const Tensor<T>& rhs, const double& thresh, const TensorType& tt) {
                        *this=gentensorT(rhs,TensorArgs(thresh,tt));
                }

                GenTensor(const SliceGenTensor<T>& rhs) : _ptr() {
                        *this=rhs;
                }

                ~GenTensor() {
                        this->clear();
                }

                gentensorT& operator=(const gentensorT& rhs) {
                        if (this != &rhs) _ptr=rhs._ptr;
                        return *this;
                }

                GenTensor& operator=(const SliceGenTensor<T>& rhs) {
                        *this=rhs._refGT.copy_slice(rhs._s);
                        return *this;
                }

                friend gentensorT copy(const gentensorT& rhs) {
                        if (rhs._ptr) return gentensorT(copy(*rhs._ptr));
                        return gentensorT();
                }


        template <class Q> operator GenTensor<Q>() const { // type conv => deep copy
            GenTensor<Q> result;
            MADNESS_EXCEPTION("no type conversion in GenTensor -- use NO_GENTENSOR",1);
            return result;
        }

                const GenTensor get_configs(const int& start, const int& end) const {
                        return gentensorT(config().get_configs(start,end));
                }

                SliceGenTensor<T> operator()(const std::vector<Slice>& s) {
                        return SliceGenTensor<T>(*this,s);
                }

                const SliceGenTensor<T> operator()(const std::vector<Slice>& s) const {
                        return SliceGenTensor<T>(*this,s);
                }

                GenTensor& operator=(const T& fac) {
                        MADNESS_EXCEPTION("no assignment with a number for GenTensor",0);
                }

        private:
                GenTensor(const Tensor<T>& t) {
                        print("in wrong constructor");
                }
        public:

                explicit GenTensor(const SRConf<T>& config) : _ptr(new configT(config)) {
                }

        private:
                gentensorT copy_slice(const std::vector<Slice>& s) const {

            // fast return if possible
            if (this->has_no_data()) {
                int k_new=s[0].end-s[0].start+1;
                return gentensorT (this->tensor_type(),k_new,this->dim());
            }

                        // consistency check
                        MADNESS_ASSERT(s.size()==this->dim());
                        MADNESS_ASSERT(s[0].step==1);

                        // fast return for full rank tensors
                        if (tensor_type()==TT_FULL) {
                                tensorT a=copy(full_tensor()(s));
                                return gentensorT(configT(a));
                        }

                        MADNESS_ASSERT(_ptr->has_structure());
//                      _ptr->make_structure();

                        // get dimensions
                        const TensorType tt=this->tensor_type();
                        const int merged_dim=this->_ptr->dim_per_vector();
                        const int dim_eff=this->_ptr->dim_eff();
                        const int rank=this->rank();
                        int k_new=s[0].end-s[0].start+1;
                        if (s[0].end<0) k_new+=this->get_k();

                        // get and reshape the vectors, slice and re-reshape again;
                        // this is shallow
                        const gentensorT& sr=*this;

                        std::vector<Tensor<T> > vectors(dim_eff,Tensor<T>());

                        for (int idim=0; idim<dim_eff; idim++) {

                                // assignment from/to slice is deep-copy
                                if (merged_dim==1) {
                                        if (rank>0) {
                                                vectors[idim]=copy(sr._ptr->ref_vector(idim)(Slice(0,rank-1),s[idim]));
                                        } else {
                                                vectors[idim]=Tensor<T>(0,s[idim].end-s[idim].start+1);
                                        }
                                } else if (merged_dim==2) {
                                        if (rank>0) {
                                                vectors[idim]=copy(sr._ptr->ref_vector(idim)(Slice(0,rank-1),s[2*idim],s[2*idim+1]));
                                        } else {
                                                vectors[idim]=tensorT(0,s[2*idim].end-s[2*idim].start+1,
                                                                                                s[2*idim+1].end-s[2*idim+1].start+1);
                                        }
                                } else if (merged_dim==3) {
                                        if (rank>0) {
                                                vectors[idim]=copy(sr._ptr->ref_vector(idim)(Slice(0,rank-1),
                                                                s[3*idim],s[3*idim+1],s[3*idim+2]));
                                        } else {
                                                vectors[idim]=tensorT(0,s[3*idim].end-s[3*idim].start+1,
                                                                s[3*idim+1].end-s[3*idim+1].start+1,
                                                                s[3*idim+2].end-s[3*idim+2].start+1);

                                        }
                                } else MADNESS_EXCEPTION("unknown number of dimensions in GenTensor::copy_slice()",0);
                        }

                        // work-around for rank==0
                        Tensor<double> weights;
                        if (rank>0) {
                                weights=copy(this->_ptr->weights_(Slice(0,rank-1)));
                        } else {
                                weights=Tensor<double>(int(0));
                        }
                        const configT conf(weights,vectors,this->dim(),k_new,tt);

                        return gentensorT(conf);

                }

        public:

                const tensorT& full_tensor() const {
                        MADNESS_ASSERT(tensor_type()==TT_FULL);
                        return _ptr->ref_vector(0);
                }

                tensorT& full_tensor() {
                        MADNESS_ASSERT(tensor_type()==TT_FULL);
                        return _ptr->ref_vector(0);
                }

                gentensorT& operator+=(const gentensorT& rhs) {

                    if (rhs.has_no_data()) return *this;
                    if (this->has_no_data()) {
                        *this=copy(rhs);
                        return *this;
                    }
                    this->gaxpy(1.0,rhs,1.0);
                        return *this;
                }

                gentensorT& operator-=(const gentensorT& rhs) {

                    if (rhs.has_no_data()) return *this;
                    if (this->has_no_data()) {
                        *this=copy(rhs);
                        return *this;
                    }
                    this->gaxpy(1.0,rhs,-1.0);
                        return *this;
                }

                gentensorT& operator+=(const SliceGenTensor<T>& rhs) {
                        const std::vector<Slice> s(this->ndim(),Slice(0,get_k()-1,1));
                        this->_ptr->inplace_add(*rhs._refGT._ptr,s,rhs._s,1.0,1.0);
                        return *this;
                }

                gentensorT& operator-=(const SliceGenTensor<T>& rhs) {
                        const std::vector<Slice> s(this->ndim(),Slice(0,get_k()-1,1));
                        this->_ptr->inplace_add(*rhs._refGT._ptr,s,rhs._s,1.0,-1.0);
                        return *this;
                }

                template<typename Q>
                GenTensor<TENSOR_RESULT_TYPE(T,Q)> operator*(const Q& x) const {
                        GenTensor<TENSOR_RESULT_TYPE(T,Q)> result(copy(*this));
                result.scale(x);
                return result;
                }

            gentensorT& gaxpy(const T alpha, const gentensorT& rhs, const T beta) {
                        MADNESS_ASSERT(this->tensor_type()==rhs.tensor_type());

                if (tensor_type()==TT_FULL) {
                full_tensor().gaxpy(alpha,rhs.full_tensor(),beta);
                return *this;
            }
                if (not (alpha==1.0)) this->scale(alpha);
                rhs.append(*this,beta);
                return *this;
            }

            template<typename Q>
            GenTensor<TENSOR_RESULT_TYPE(T,Q)>& scale(const Q& dfac) {
                        if (!_ptr) return *this;
                        if (tensor_type()==TT_FULL) {
                                full_tensor().scale(dfac);
                        } else {
                                _ptr->scale(dfac);
                        }
                        return *this;
                };

                void normalize() {
                        if (tensor_type()==TT_2D) _ptr->normalize();
                };

                void fillrandom(const int r=1) {
                        if (tensor_type()==TT_FULL) full_tensor().fillrandom();
                        else _ptr->fillWithRandom(r);
                }

                bool has_data() const {
                    if (_ptr) return true;
                    return false;
                }

                bool has_no_data() const {return (!has_data());}

                long rank() const {
                        if (_ptr) return _ptr->rank();
                        else return 0;
                };

                unsigned int dim() const {return _ptr->dim();};

                long dim(const int& i) const {return _ptr->get_k();};

                long ndim() const {
                        if (_ptr) return _ptr->dim();
                        return -1;
                };

                unsigned int get_k() const {return _ptr->get_k();};

                TensorType tensor_type() const {
                        if (_ptr) return _ptr->type();
                        return TT_NONE;
                };

        std::string what_am_i() const {return TensorArgs::what_am_i(_ptr->type());};

                size_t size() const {
                        if (_ptr) return _ptr->nCoeff();
                        return 0;
                };

                size_t real_size() const {
                        if (_ptr) return _ptr->real_size()+sizeof(*this);
                        return 0;
                };

                double normf() const {
                        if (has_no_data()) return 0.0;
                        return config().normf();
                };

        double svd_normf() const {
            if (has_no_data()) return 0.0;
            if (tensor_type()==TT_2D) return config().svd_normf();
            return config().normf();
        };

                T trace(const GenTensor<T>& rhs) const {
                        return this->trace_conj(rhs);
                }

                template<typename Q>
                TENSOR_RESULT_TYPE(T,Q) trace_conj(const GenTensor<Q>& rhs) const {
            if (TensorTypeData<T>::iscomplex) MADNESS_EXCEPTION("no complex trace in GenTensor, sorry",1);
            if (TensorTypeData<Q>::iscomplex) MADNESS_EXCEPTION("no complex trace in GenTensor, sorry",1);

            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
                        // fast return if possible
                        if ((this->rank()==0) or (rhs.rank()==0)) return resultT(0.0);

                        MADNESS_ASSERT(compatible(*this,rhs));
                        MADNESS_ASSERT(this->tensor_type()==rhs.tensor_type());

                        return overlap(*(this->_ptr),*rhs._ptr);
                }

                template<typename Q>
                TENSOR_RESULT_TYPE(T,Q) trace_conj(const Tensor<Q>& rhs) const {
            if (TensorTypeData<T>::iscomplex) MADNESS_EXCEPTION("no complex trace in GenTensor, sorry",1);
            if (TensorTypeData<Q>::iscomplex) MADNESS_EXCEPTION("no complex trace in GenTensor, sorry",1);

            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
                        // fast return if possible
                        if ((this->rank()==0)) return resultT(0.0);

                        // fast return if this is a full tensor
                        // otherwise reconstruct this to a full tensor, since it's presumably
                        // faster than turning the full tensor rhs into a low-rank approximation
                        if (this->tensor_type()==TT_FULL) return full_tensor().trace_conj(rhs);
                        else return full_tensor_copy().trace_conj(rhs);

                }

        GenTensor<T>& emul(const GenTensor<T>& t) {
                print("no GenTensor<T>::emul yet");
                MADNESS_ASSERT(0);
        }

                Tensor<T> full_tensor_copy() const {
                        const TensorType tt=tensor_type();
                        if (tt==TT_NONE) return Tensor<T>();
                        else if (tt==TT_2D) return this->reconstruct_tensor();
                        else if (tt==TT_FULL) {
                                return copy(full_tensor());
                        } else {
                                print(TensorArgs::what_am_i(tt),"unknown tensor type");
                                MADNESS_ASSERT(0);
                        }
                }

                void reduce_rank(const double& eps) {

                        if (rank()==0) return;

                        // direct reduction on the polynomial values on the Tensor
                        if (tensor_type()==TT_FULL or tensor_type()==TT_NONE) {
                                return;
                        } else if (this->tensor_type()==TT_2D) {
                                config().divide_and_conquer_reduce(eps*facReduce());
                        } else {
                                MADNESS_EXCEPTION("unknown tensor type in GenTensor::reduceRank()",0);
                        }
                        MADNESS_ASSERT(this->_ptr->has_structure() or this->rank()==0);
                }

                void printCoeff(const std::string title) const {
                        print("printing SepRep",title);
                        print(_ptr->weights_);
                        for (unsigned int idim=0; idim<this->_ptr->dim_eff(); idim++) {
                                print("coefficients for dimension",idim);
                                print(_ptr->vector_[idim]);
                        }
                }

                Tensor<T> reconstruct_tensor() const {

                        /*
                         * reconstruct the tensor first to the configurational dimension,
                         * then to the real dimension
                         */

                        // fast return for full rank tensors
                        if (tensor_type()==TT_FULL) return full_tensor();

                        // for convenience
                        const unsigned int conf_dim=this->_ptr->dim_eff();
                        const unsigned int conf_k=this->_ptr->kVec();                   // possibly k,k*k,..
                        const long rank=this->rank();
                        long d[TENSOR_MAXDIM];

                        // fast return if possible
                        if (rank==0) {
                                // reshape the tensor to the really required one
                                const long k=this->get_k();
                                const long dim=this->dim();
                                for (long i=0; i<dim; i++) d[i] = k;

                                return Tensor<T> (dim,d,true);
                        }


                        // set up result Tensor (in configurational dimensions)
                        for (long i=0; i<conf_dim; i++) d[i] = conf_k;
                        tensorT s(conf_dim,d,true);

                        // flatten this
                        gentensorT sr=*this;

                        // and a scratch Tensor
                        Tensor<T>  scr(rank);
                        Tensor<T>  scr1(rank);
                        Tensor<T>  scr2(rank);

                        if (conf_dim==1) {

                                for (unsigned int i0=0; i0<conf_k; i0++) {
                                        scr=sr._ptr->weights_;
                                        //                                      scr.emul(F[0][i0]);
                                        T buffer=scr.sum();
                                        s(i0)=buffer;
                                }

                        } else if (conf_dim==2) {


                                //                              tensorT weight_matrix(rank,rank);
                                //                              for (unsigned int r=0; r<rank; r++) {
                                //                                      weight_matrix(r,r)=this->weight(r);
                                //                              }
                                //                              s=inner(weight_matrix,sr._ptr->refVector(0));
                                //                              s=inner(s,sr._ptr->refVector(1),0,0);
//                tensorT sscr=copy(sr._ptr->ref_vector(0)(sr._ptr->c0()));
                                tensorT sscr=copy(sr._ptr->flat_vector(0));
                                for (unsigned int r=0; r<rank; r++) {
                                        const double w=_ptr->weights(r);
                                        for (unsigned int k=0; k<conf_k; k++) {
                                                sscr(r,k)*=w;
                                        }
                                }
                                inner_result(sscr,sr._ptr->flat_vector(1),0,0,s);


                        } else if (conf_dim==3) {

                MADNESS_EXCEPTION("flat/structure issue",1);
                                for (unsigned int i0=0; i0<conf_k; i0++) {
                                        scr=copy(sr._ptr->weights_(Slice(0,sr.rank()-1)));
                                        scr.emul(sr._ptr->ref_vector(0)(Slice(0,rank-1),i0));
                                        for (unsigned int i1=0; i1<conf_k; i1++) {
                                                scr1=copy(scr);
                                                scr1.emul(sr._ptr->ref_vector(1)(Slice(0,rank-1),i1));
                                                for (unsigned int i2=0; i2<conf_k; i2++) {
                                                        scr2=copy(scr1);
                                                        scr2.emul(sr._ptr->ref_vector(2)(Slice(0,rank-1),i2));
                                                        s(i0,i1,i2)=scr2.sum();
                                                }
                                        }
                                }
                        } else {
                                print("only config_dim=1,2,3 in GenTensor::reconstructTensor");
                                MADNESS_ASSERT(0);
                        }


                        // reshape the tensor to the really required one
                        const long k=this->get_k();
                        const long dim=this->dim();
                        for (long i=0; i<dim; i++) d[i] = k;

                        Tensor<T> s2=s.reshape(dim,d);
                        return s2;
                }

                void append(gentensorT& rhs, const T fac=1.0) const {
                        rhs.config().append(*this->_ptr,fac);
                }

                void add_SVD(const gentensorT& rhs, const double& thresh) {
                        if (rhs.has_no_data()) return;
            if (has_no_data()) {
                *this=rhs;
                return;
            }
                        if (tensor_type()==TT_FULL or tensor_type()==TT_NONE) {
                                this->full_tensor()+=rhs.full_tensor();
                                return;
                        }
                        config().add_SVD(rhs.config(),thresh*facReduce());
                }

                friend bool compatible(const gentensorT& rhs, const gentensorT& lhs) {
                        return ((rhs.tensor_type()==lhs.tensor_type()) and (rhs.get_k()==lhs.get_k())
                                        and (rhs.dim()==lhs.dim()));
                };

                gentensorT transform(const Tensor<T> c) const {
//                      _ptr->make_structure();
                    if (has_no_data()) return gentensorT();
                        MADNESS_ASSERT(_ptr->has_structure());
                        return gentensorT (this->_ptr->transform(c));
                }

                template<typename Q>
                gentensorT general_transform(const Tensor<Q> c[]) const {
//                  this->_ptr->make_structure();
                    if (has_no_data()) return gentensorT();
            MADNESS_ASSERT(_ptr->has_structure());
                        return gentensorT (this->config().general_transform(c));
                }

                gentensorT transform_dir(const Tensor<T>& c, const int& axis) const {
//            this->_ptr->make_structure();
            MADNESS_ASSERT(_ptr->has_structure());
            return GenTensor<T>(this->_ptr->transform_dir(c,axis));
                }

                const SRConf<T>& config() const {return *_ptr;}

                SRConf<T>& config() {return *_ptr;}

                static double fac_reduce() {return facReduce();};

        private:

                void clear() {_ptr.reset();};

                void inplace_add(const gentensorT& rhs, const std::vector<Slice>& lhs_s,
                                const std::vector<Slice>& rhs_s, const double alpha, const double beta) {

                        // fast return if possible
                        if (rhs.has_no_data() or rhs.rank()==0) return;

                        if (this->has_data()) MADNESS_ASSERT(this->tensor_type()==rhs.tensor_type());

                        // no fast return possible!!!
                        //                      if (this->rank()==0) {
                        //                              // this is a deep copy
                        //                              *this=rhs(rhs_s);
                        //                              this->scale(beta);
                        //                              return;
                        //                      }

                        if (tensor_type()==TT_FULL) {
                                full_tensor()(lhs_s).gaxpy(alpha,rhs.full_tensor()(rhs_s),beta);
                        } else {
//                              rhs._ptr->make_structure();
//                              _ptr->make_structure();
                                MADNESS_ASSERT(_ptr->has_structure());
                MADNESS_ASSERT(rhs._ptr->has_structure());
                                this->_ptr->inplace_add(*rhs._ptr,lhs_s,rhs_s, alpha, beta);
                        }
                }

                void computeSVD(const double& eps,const Tensor<T>& values_eff) {

                        // SVD works only with matrices (2D)
                        MADNESS_ASSERT(values_eff.ndim()==2);
                        MADNESS_ASSERT(this->tensor_type()==TT_2D);

                        // fast return if possible
                        if (values_eff.normf()<eps*facReduce()) {
                                _ptr=sr_ptr(new configT(_ptr->dim(),_ptr->get_k(),tensor_type()));
                                return;
                        }

                        // output from svd
                        Tensor<T> U;
                        Tensor<T> VT;
                        Tensor< typename Tensor<T>::scalar_type > s;

                        svd(values_eff,U,s,VT);

                        // find the maximal singular value that's supposed to contribute
                        // singular values are ordered (largest first)
                        const double thresh=eps*facReduce();
                        long i=SRConf<T>::max_sigma(thresh,s.dim(0),s);

                        // convert SVD output to our convention
                        if (i>=0) {
                            // copy to have contiguous and tailored singular vectors
                                _ptr=sr_ptr(new configT(copy(s(Slice(0,i))), copy(transpose(U(_,Slice(0,i)))),
                                        copy(VT(Slice(0,i),_)), dim(), get_k()));
                MADNESS_ASSERT(this->_ptr->get_k()==this->_ptr->vector_[0].dim(1));
                MADNESS_ASSERT(this->_ptr->rank()==this->_ptr->vector_[0].dim(0));
                MADNESS_ASSERT(this->_ptr->rank()==this->_ptr->weights_.dim(0));

                        } else {
                                _ptr=sr_ptr(new configT(dim(),get_k(),tensor_type()));
                        }
                }
        };



        template <typename T>
        class SliceGenTensor {

        private:
                friend class GenTensor<T>;

                GenTensor<T>& _refGT;
                std::vector<Slice> _s;

                // all ctors are private, only accessible by GenTensor

                SliceGenTensor<T> () {}

                SliceGenTensor<T> (const GenTensor<T>& gt, const std::vector<Slice>& s)
                                : _refGT(const_cast<GenTensor<T>& >(gt))
                                , _s(s) {}

        public:

                SliceGenTensor<T>& operator=(const GenTensor<T>& rhs) {
                        print("You don't want to assign to a SliceGenTensor; use operator+= instead");
                        MADNESS_ASSERT(0);
                        return *this;
                };

                SliceGenTensor<T>& operator=(const SliceGenTensor<T>& rhs) {
                        print("You don't want to assign to a SliceGenTensor; use operator+= instead");
                        MADNESS_ASSERT(0);
                        return *this;
                };

                SliceGenTensor<T>& operator+=(const GenTensor<T>& rhs) {
                        std::vector<Slice> s(this->_refGT.ndim(),Slice(_));
                        _refGT.inplace_add(rhs,this->_s,s,1.0,1.0);
                        return *this;
                }

                SliceGenTensor<T>& operator-=(const GenTensor<T>& rhs) {
                        std::vector<Slice> s(this->_refGT.ndim(),Slice(_));
                        _refGT.inplace_add(rhs,this->_s,s,1.0,-1.0);
                        return *this;
                }

                SliceGenTensor<T>& operator+=(const SliceGenTensor<T>& rhs) {
                        _refGT.inplace_add(GenTensor<T>(*rhs._refGT._ptr),this->_s,rhs._s,1.0,1.0);
                        return *this;
                }

                SliceGenTensor<T>& operator=(const T& number) {
                        MADNESS_ASSERT(number==T(0.0));
                        if (this->_refGT.tensor_type()==TT_FULL) {
                                _refGT.full_tensor()(_s)=0.0;
                        } else {
                                const GenTensor<T>& tmp=(this->_refGT);
                                _refGT.inplace_add(tmp,_s,_s,1.0,-1.0);
                        }
                        return *this;
                }

                friend GenTensor<T> copy(const SliceGenTensor<T>& rhs) {
                    if (rhs._refGT.has_data()) return GenTensor<T>(rhs);
                    return GenTensor<T>();
                }

        };


    template <class T, class Q>
    GenTensor< TENSOR_RESULT_TYPE(T,Q) > transform(const GenTensor<Q>& t, const Tensor<T>& c) {
        return t.transform(c);
    }

    struct RankReduceWrapper {
        double eps;
        RankReduceWrapper(const double& e) : eps(e) {}

        template<typename T>
        void operator()(GenTensor<T>& g) const {
                g.reduce_rank(eps);
        }
    };

    template<typename T>
    bool compare_rank(const GenTensor<T>& g1, const GenTensor<T>& g2) {
        return g1.rank()<g2.rank();
    }



    template<typename T>
        GenTensor<T> reduce(std::list<GenTensor<T> >& addends, double eps, bool are_optimal=false) {

        // fast return
        if (addends.size()==0) return GenTensor<T>();

        typedef typename std::list<GenTensor<T> >::iterator iterT;

        // make error relative
        eps=eps/addends.size();

        // if the addends are not in SVD format do that now so that we can call add_svd later
        if (not are_optimal) {
                RankReduceWrapper f(eps);
                std::for_each(addends.begin(),addends.end(),f);
        }

        // remove zero ranks and sort the list according to the gentensor's ranks
                addends.remove_if(has_zero_rank<T>);
        if (addends.size()==0) return GenTensor<T>();
                addends.sort(compare_rank<T>);

        // do the additions
        GenTensor<T> result=copy(addends.front());
        for (iterT it=++addends.begin(); it!=addends.end(); ++it) {
                GenTensor<T>& rhs=*it;
                MADNESS_ASSERT(&rhs!=&result);
                result.add_SVD(rhs,eps);
        }
        addends.clear();

        return result;
    }


    template <class T, class Q>
    GenTensor<TENSOR_RESULT_TYPE(T,Q)> transform_dir(const GenTensor<Q>& t, const Tensor<T>& c,
                                          int axis) {

        return t.transform_dir(c,axis);
    }

    template<typename T>
    static inline
    std::ostream& operator<<(std::ostream& s, const GenTensor<T>& g) {
        std::string str="GenTensor has no data";
        if (g.has_no_data()) {
            std::string str="GenTensor has no data";
            s << str.c_str() ;
        } else {
            str="GenTensor has data";
            s << str.c_str() << g.config();
        }
        return s;
    }

    namespace archive {
        template <class Archive, typename T>
                struct ArchiveStoreImpl< Archive, GenTensor<T> > {

                        friend class GenTensor<T>;
                        static void store(const Archive& ar, const GenTensor<T>& t) {
                                bool exist=t.has_data();
                                ar & exist;
                                if (exist) ar & t.config();
                        };
                };


                template <class Archive, typename T>
                struct ArchiveLoadImpl< Archive, GenTensor<T> > {

                        friend class GenTensor<T>;
                        static void load(const Archive& ar, GenTensor<T>& t) {
                                // check for pointer existence
                                bool exist=false;
                                ar & exist;
                                if (exist) {
                                        SRConf<T> conf;
                                        ar & conf;
                                        //t.config()=conf;
                                        t=GenTensor<T>(conf);
                                }
                        };
                };
    };

     template <class T, class Q>
     GenTensor<TENSOR_RESULT_TYPE(T,Q)> outer(const GenTensor<T>& lhs2, const GenTensor<Q>& rhs2) {
        return outer_low_rank(lhs2.full_tensor(),rhs2.full_tensor());
     }

     template <class T, class Q>
     GenTensor<TENSOR_RESULT_TYPE(T,Q)> outer_low_rank(const Tensor<T>& lhs2, const Tensor<Q>& rhs2) {

        typedef TENSOR_RESULT_TYPE(T,Q) resultT;

        // srconf is shallow, do deep copy here
        const Tensor<T> lhs=copy(lhs2);
        const Tensor<Q> rhs=copy(rhs2);

        const long k=lhs.dim(0);
        const long dim=lhs.ndim()+rhs.ndim();
        long size=1;
        for (int i=0; i<lhs.ndim(); ++i) size*=k;
        MADNESS_ASSERT(size==lhs.size());
        MADNESS_ASSERT(size==rhs.size());
        MADNESS_ASSERT(lhs.size()==rhs.size());

                Tensor<double> weights(1);
                weights=1.0;

                SRConf<resultT> srconf(weights,lhs.reshape(1,lhs.size()),rhs.reshape(1,rhs.size()),dim,k);
                GenTensor<resultT> coeff(srconf);
                coeff.normalize();
                return coeff;
     }

    #endif /* HAVE_GENTENSOR */

    template<typename T>
    void change_tensor_type(GenTensor<T>& t, const TensorArgs& targs) {

        // fast return if possible
        const TensorType current_type=t.tensor_type();
        if (current_type==targs.tt) return;
        if (t.has_no_data()) return;

        // for now
        MADNESS_ASSERT(targs.tt==TT_FULL or targs.tt==TT_2D);
        MADNESS_ASSERT(current_type==TT_FULL or current_type==TT_2D);

        GenTensor<T> result;
        if (targs.tt==TT_FULL) {
            result=GenTensor<T>(t.full_tensor_copy(),targs);
        } else if (targs.tt==TT_2D) {
            MADNESS_ASSERT(current_type==TT_FULL);
            result=GenTensor<T>(t.full_tensor(),targs);
        }

        t=result;

    }


    template <class T, class Q>
    GenTensor<TENSOR_RESULT_TYPE(T,Q)> general_transform(const GenTensor<T>& t, const Tensor<Q> c[]) {
        return t.general_transform(c);
    }

    template <class T>
    GenTensor<T> general_transform(const GenTensor<T>& t, const Tensor<T> c[]) {
        return t.general_transform(c);
    }

    template <typename T, typename Q>
    typename IsSupported < TensorTypeData<Q>, GenTensor<T> >::type
    operator*(const Q& x, const GenTensor<T>& t) {
        return t*x;
    }



}

#endif /* GENTENSOR_H_ */