Commit e5fefdc7 authored by Alberto's avatar Alberto
Browse files

LinearPrediction: fixed

parent 7f703dc5
......@@ -30,7 +30,7 @@ typedef PsimagLite::LinearPrediction<FieldType> LinearPredictionType;
void usage(const char *progName)
{
std::cerr<<"Usage: "<<progName<<" -f file -l label -p n\n";
std::cerr<<"Usage: "<<progName<<" -f file -l label -p n -q q[>1]\n";
}
int main(int argc,char *argv[])
......@@ -39,8 +39,9 @@ int main(int argc,char *argv[])
PsimagLite::String file="";
PsimagLite::String label="";
SizeType p = 0;
SizeType q = 2;
while ((opt = getopt(argc, argv,
"f:p:l:")) != -1) {
"f:p:l:q:")) != -1) {
switch (opt) {
case 'f':
file = optarg;
......@@ -51,13 +52,16 @@ int main(int argc,char *argv[])
case 'p':
p = atoi(optarg);
break;
case 'q':
q = atoi(optarg);
break;
default:
usage(argv[0]);
return 1;
}
}
// sanity checks:
if (file=="" || label=="" || p==0) {
if (file=="" || label=="" || p==0 || q<2) {
usage(argv[0]);
return 1;
}
......@@ -67,12 +71,19 @@ int main(int argc,char *argv[])
io.read(y,label);
SizeType n = y.size();
std::cout<<"#Found "<<n<<" points in file "<<file<<"\n";
LinearPredictionType linearPrediction(y);
linearPrediction.predict(p);
LinearPredictionType linearPrediction(y,q);
linearPrediction.predict(q);
for (SizeType i=0;i<p;i++) {
linearPrediction.linearPredictionfunction(y,q);
linearPrediction.predict(q);
}
for (SizeType i=0;i<p+n;i++) {
std::cout<<i<<" "<<linearPrediction(i)<<"\n";
std::cout<<linearPrediction(i)<<"\n";
}
}
// Test Function in series.txt is: f[i] = 0.5*(cos((i-1)*pi/100)+cos((i-1)*pi/20))*exp(-(i-1)/100)
#Series
20
Series
50
1
1.0024222
1.0088111
1.0167278
1.0228333
1.0237833
1.0171778
1.0023944
0.98105
0.95696112
0.93562778
0.92321112
0.92525
0.94531666
0.98386666
1.0375833
1.0993333
1.1588056
1.20385
1.2222778
1.2038722
1.1422722
1.0364556
0.89142778
0.71813334
0.53235778
0.35288222
0.19909111
0.088361112
0.033652778
0.041651166
0.11176167
0.23611166
0.40057666
0.58667778
0.77406112
0.94321666
1.0780056
1.1676111
1.2077333
1.2007167
1.1548278
1.0826056
0.99881666
0.91811666
0.85295556
0.81197778
0.79910556
0.81345556
0.85000556
0.9837109906618
0.9552444053212
0.9154059818405
0.8652541558705
0.8060694941695
0.7393192992506
0.6666183469137
0.5896867938457
0.5103063412053
0.4302757613054
0.3513668881767
0.2752821396547
0.2036145800083
0.1378114499733
0.07914198783209
0.02867024385213
-0.01276654573324
-0.04457360629771
-0.06640802967052
-0.07818183001025
-0.08005870195437
-0.07244468542249
-0.05597308654312
-0.03148413783927
-4.32321280481e-17
0.03730419034575
0.07913243961241
0.1241023338419
0.1707809038286
0.2177210123731
0.2634973624311
0.3067412520585
0.3461732500491
0.3806330334234
0.4091057124367
0.4307440681104
0.4448862387402
0.4510685124194
0.4490330092085
0.4387301659654
0.4203160658007
0.394144779489
0.3607560049616
0.3208584004382
0.2753091043404
0.225090018715
0.1712815007022
0.1150341572763
0.0575394711582
......@@ -39,7 +39,7 @@ must include the following acknowledgment:
"This product includes software produced by UT-Battelle,
LLC under Contract No. DE-AC05-00OR22725 with the
Department of Energy."
*********************************************************
DISCLAIMER
......@@ -79,100 +79,176 @@ DISCLOSED WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.
* Extrapolating a "time" series
* see extrapolation.tex for more details
*/
#ifndef LINEAR_PREDICTION_H
#define LINEAR_PREDICTION_H
#include "Matrix.h"
#include "LAPACK.h"
#include "BLAS.h"
#include <complex>
#ifdef USE_GSL
extern "C" {
#include <gsl/gsl_poly.h>
}
#endif
namespace PsimagLite {
template<typename FieldType>
class LinearPrediction {
typedef Matrix<FieldType> MatrixType;
public:
LinearPrediction(const typename Vector<FieldType>::Type& y)
: y_(y)
{
SizeType ysize = y.size();
if (ysize&1) throw RuntimeError(
"LinearPrediction::ctor(...): data set must contain an even number of points\n");
SizeType n = ysize/2;
MatrixType A(n,n);
typename Vector<FieldType>::Type B(n);
computeA(A);
computeB(B);
computeD(A,B);
template<typename FieldType>
class LinearPrediction {
typedef Matrix<FieldType> MatrixType;
public:
LinearPrediction(const typename Vector<FieldType>::Type& y, SizeType p)
: y_(y), p_(p)
{
SizeType ysize = y.size();
if (ysize&1) throw RuntimeError(
"LinearPrediction::ctor(...): data set must contain an even number of points\n");
SizeType n = ysize/2;
MatrixType A(p,p);
typename Vector<FieldType>::Type B(p);
computeA(A,n);
computeB(B,n);
computeD(A,B);
}
const FieldType& operator()(SizeType i) const
{
return y_[i];
}
void linearPredictionfunction(const typename Vector<FieldType>::Type& y, SizeType p)
{
SizeType ysize = y.size();
if (ysize&1) throw RuntimeError(
"LinearPrediction::ctor(...): data set must contain an even number of points\n");
SizeType n = ysize/2;
MatrixType A(p,p);
typename Vector<FieldType>::Type B(p);
computeA(A,n);
computeB(B,n);
computeD(A,B);
}
void predict(SizeType p)
{
SizeType n = y_.size();
for (SizeType i=n;i<n+p;i++) {
FieldType sum = 0;
for (SizeType j=0;j<d_.size();j++) {
sum += d_[j] * y_[i-j-1];
}
y_.push_back(sum);
}
}
void predict2(SizeType p)
{
// Fix roots
typedef std::complex<double> Complex;
SizeType twicep = 2*p;
SizeType pp1 = p+1;
typename Vector<FieldType>::Type nd(p),aa(pp1),zz(twicep);
std::vector<Complex> roots(p);
std::vector<Complex> ab(pp1);
const FieldType& operator()(SizeType i) const
{
return y_[i];
aa[p]= 1;
for (SizeType j=0;j<p;j++) {
aa[j]= -d_[p-1-j];
}
void predict(SizeType p)
{
SizeType n = y_.size();
for (SizeType i=n;i<n+p;i++) {
FieldType sum = 0;
for (SizeType j=0;j<d_.size();j++)
sum += d_[j]*y_[i-j-1];
y_.push_back(sum);
}
gsl_poly_complex_workspace* w = gsl_poly_complex_workspace_alloc(pp1);
gsl_poly_complex_solve (&aa[0], pp1, w, &zz[0]);
gsl_poly_complex_workspace_free(w);
for (SizeType j=0;j<p;j++) {
roots[j] = Complex(zz[2*j],zz[2*j+1]);
}
private:
//! Note: A and B cannot be const. here due to the ultimate
//! call to BLAS::GEMV
void computeD(MatrixType& A,typename Vector<FieldType>::Type& B)
{
SizeType n = B.size();
typename Vector<int>::Type ipiv(n); // use signed integers here!!
int info = 0;
psimag::LAPACK::GETRF(n, n, &(A(0,0)), n, &(ipiv[0]), info);
typename Vector<FieldType>::Type work(2);
int lwork = -1; // query mode
psimag::LAPACK::GETRI(n, &(A(0,0)), n, &(ipiv[0]),
&(work[0]), lwork,info );
lwork = static_cast<int>(work[0]);
if (lwork<=0) throw
RuntimeError("LinearPrediction:: internal error\n");
work.resize(lwork);
// actual work:
psimag::LAPACK::GETRI(n, &(A(0,0)), n, &(ipiv[0]),
&(work[0]), lwork,info );
d_.resize(n);
psimag::BLAS::GEMV('N',n,n,1.0,&(A(0,0)),n,&(B[0]),1,0.0,&(d_[0]),1);
for (SizeType j=0;j<p;j++) {
//Look for a root outside the unit circle, and put it to 0
if (abs(roots[j]) > 1.0) {
roots[j]= 0.0;//roots[j]/abs(roots[j]);
}
}
//Now reconstruct the polynomial coefficients
ab[0] = -roots[0];
ab[1] = 1.0;
for (SizeType j=1;j<p;j++) {
ab[j+1] = 1.0;
for (SizeType i=j;i>=1;i--) {
ab[i]=ab[i-1]-roots[j]*ab[i];
}
ab[0]= -roots[j]*ab[0];
}
for (SizeType j=0;j<p;j++) {
nd[p-1-j] = -PsimagLite::real(ab[j]);
}
void computeA(MatrixType& A) const
{
SizeType n = A.rows();
for (SizeType l=0;l<n;l++) {
for (SizeType j=0;j<n;j++) {
A(l,j) = 0;
for (SizeType i=n;i<2*n;i++)
A(l,j) += y_[i-l-1] * y_[i-j-1];
}
SizeType n = y_.size();
y_.resize(n+p);
for (SizeType i=n;i<n+p;i++) {
FieldType sum = 0;
for (SizeType j=0;j<d_.size();j++) {
sum += nd[j] * y_[i-j-1];
}
y_[i] = sum;
}
}
private:
//! Note: A and B cannot be const. here due to the ultimate
//! call to BLAS::GEMV
void computeD(MatrixType& A,typename Vector<FieldType>::Type& B)
{
SizeType p = B.size();
typename Vector<int>::Type ipiv(p); // use signed integers here!!
int info = 0;
psimag::LAPACK::GETRF(p, p, &(A(0,0)), p, &(ipiv[0]), info);
void computeB(typename Vector<FieldType>::Type& B) const
{
SizeType n = B.size();
for (SizeType l=0;l<n;l++) {
B[l] = 0;
typename Vector<FieldType>::Type work(2);
int lwork = -1; // query mode
psimag::LAPACK::GETRI(p, &(A(0,0)), p, &(ipiv[0]),
&(work[0]), lwork,info );
lwork = static_cast<int>(work[0]);
if (lwork<=0) throw
RuntimeError("LinearPrediction:: internal error\n");
work.resize(lwork);
// actual work:
psimag::LAPACK::GETRI(p, &(A(0,0)), p, &(ipiv[0]),
&(work[0]), lwork,info );
d_.resize(p);
psimag::BLAS::GEMV('N',p,p,1.0,&(A(0,0)),p,&(B[0]),1,0.0,&(d_[0]),1);
}
void computeA(MatrixType& A, SizeType n) const
{
SizeType p = A.rows();
for (SizeType l=0;l<p;l++) {
for (SizeType j=0;j<p;j++) {
A(l,j) = 0;
for (SizeType i=n;i<2*n;i++)
B[l] += y_[i-l-1] * y_[i];
A(l,j) += y_[i-l-1] * y_[i-j-1];
}
}
}
void computeB(typename Vector<FieldType>::Type& B, SizeType n) const
{
SizeType p = B.size();
for (SizeType l=0;l<p;l++) {
B[l] = 0;
for (SizeType i=n;i<2*n;i++)
B[l] += y_[i-l-1] * y_[i];
}
}
typename Vector<FieldType>::Type y_;
typename Vector<FieldType>::Type d_;
}; // class LinearPrediction
typename Vector<FieldType>::Type y_;
SizeType p_;
typename Vector<FieldType>::Type d_;
}; // class LinearPrediction
} // namespace PsimagLite
/*@}*/
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment