spectrum.cc
//
// Distributed under the ITensor Library License, Version 1.2
// (See accompanying LICENSE file.)
//
#include <algorithm>
#include <utility>
#include "itensor/spectrum.h"
using std::pair;
using std::make_pair;
namespace itensor {
struct OrderSecond
{
template<typename PairType>
bool
operator()(PairType const& i, PairType const& j) const
{
return i.second > j.second;
}
};
Spectrum::
Spectrum(Args const& args)
:
truncerr_(NAN)
{
computeTruncerr(args);
}
//Spectrum::
//Spectrum(ITensor const& D, Args const& args)
// :
// truncerr_(0)
// {
// Error("Spectrum ITensor constructor not yet implemented");
// //if(D.type() != ITensor::Diag) Error("Spectrum may only be constructed from Diag type ITensor.");
// //eigs_ = D.diag();
// //for(auto n = 1l; n <= eigs_.size(); ++n)
// // eigs_(n) = sqr(eigs_(n));
// //computeTruncerr(args);
// }
//
//Spectrum::
//Spectrum(const IQTensor& D, const Args& args)
// :
// truncerr_(0)
// {
// Error("Spectrum IQTensor constructor not yet implemented");
// //std::vector<OrderSecond::value_type> eigs;
// //eigs.reserve(D.indices().front().m());
//
// //for(const ITensor& t : D.blocks())
// // {
// // if(t.type() != ITensor::Diag)
// // Error("Spectrum may only be constructed from IQTensor containing only Diag type ITensor.");
// // const Vector svals = t.diag();
// // const QN q = itensor::qn(D,t.indices().front());
// // for(int n = 1; n <= svals.size(); ++n)
// // {
// // eigs.push_back(std::make_pair(q,sqr(svals(n))));
// // }
// // }
// //std::sort(eigs.begin(),eigs.end(),OrderSecond());
//
// //qns_.resize(eigs.size());
// //eigs_.ReDimension(eigs.size());
// //for(size_t j = 0; j < eigs.size(); ++j)
// // {
// // qns_.at(j) = eigs.at(j).first;
// // eigs_[j] = eigs.at(j).second;
// // }
// //computeTruncerr(args);
// }
Spectrum::
Spectrum(Vector const& eigs, Args const& args)
:
eigs_(eigs)
{
computeTruncerr(args);
}
Spectrum::
Spectrum(Vector const& eigs,
QNStorage const& qns,
Args const& args)
:
eigs_(eigs),
qns_(qns)
{
computeTruncerr(args);
}
QN Spectrum::
qn(int n) const
{
if(qns_.empty()) return QN();
return qns_.at(n-1);
}
void Spectrum::
read(std::istream& s)
{
Error("Spectrum::read not currently implemented");
//s.read((char*)&truncerr_,sizeof(truncerr_));
//eigs_.read(s);
//size_t sz = 0;
//s.read((char*)&sz,sizeof(sz));
//qns_.resize(sz);
//for(size_t j = 0; j < sz; ++j)
// {
// s.read((char*)&qns_[j],sizeof(qns_[j]));
// }
}
void Spectrum::
write(std::ostream& s) const
{
Error("Spectrum::read not currently implemented");
//s.write((char*)&truncerr_,sizeof(truncerr_));
//eigs_.write(s);
//size_t sz = qns_.size();
//s.write((char*)&sz,sizeof(sz));
//for(size_t j = 0; j < sz; ++j)
// {
// s.write((char*)&qns_[j],sizeof(qns_[j]));
// }
}
void Spectrum::
computeTruncerr(Args const& args)
{
if(args.defined("Truncerr"))
{
truncerr_ = args.getReal("Truncerr");
return;
}
if(eigs_.size() > 0)
{
// Note: This only makes sense if the spectrum was normalized before truncation.
// Otherwise this will mostly return zero.
truncerr_ = 1.-sumels(eigs_);
if(truncerr_ < 0) truncerr_ = 0;
}
}
std::ostream&
operator<<(std::ostream & s, Spectrum const& spec)
{
auto& eigs = spec.eigsKept();
auto N = eigs.size();
if(N > 0)
{
decltype(N) max_show = 20;
auto stop = std::min(N,max_show);
s << " Eigs kept:";
for(auto j : range(stop))
{
s << format(eigs(j) > 1E-3 ? (" %.3f") : (" %.3E"), eigs(j));
}
s << format("\n Trunc. error = %.3E\n", spec.truncerr());
}
return s;
}
} //namespace itensor