ADMB Documentation  11.1.2503
 All Classes Files Functions Variables Typedefs Friends Defines
orthpoly.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: orthpoly.cpp 1674 2014-02-25 19:04:59Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <fvar.hpp>
00012 
00017 dmatrix orthpoly(int n,int deg)
00018 {
00019   int j; int is; int ik;
00020   dmatrix ocoff(0,deg,1,n);
00021   double sum;
00022   ocoff(0)=sqrt(double(n));
00023   for (is=1; is<=deg; is++)
00024   {
00025     for (j=1; j<=n; j++)
00026     {
00027       ocoff(is,j)=pow(double(j),is);
00028     }
00029   }
00030   for (is=0; is<=deg; is++) /* L1000  */
00031   {
00032     for (ik=0; ik<=is-1; ik++) /* L2000  */
00033     {
00034       sum=ocoff(is)*ocoff(ik);
00035       ocoff(is)-=sum*ocoff(ik);
00036     }
00037     sum=norm2(ocoff(is));
00038     ocoff(is)=ocoff(is)/sqrt(sum);
00039   }
00040   return trans(ocoff);
00041 }
00042 
00047 dmatrix orthpoly(int n,int deg,int skip)
00048 {
00049   int j; int is; int ik;
00050   dmatrix ocoff(0,deg,1,n);
00051   double sum;
00052   ocoff(0)=sqrt(double(n));
00053   for (is=1; is<=deg; is++)
00054   {
00055     for (j=1; j<=n; j++)
00056     {
00057       ocoff(is,j)=pow(double(j),is);
00058     }
00059   }
00060   for (is=0; is<=deg; is++) /* L1000  */
00061   {
00062     for (ik=0; ik<=is-1; ik++) /* L2000  */
00063     {
00064       sum=ocoff(is)*ocoff(ik);
00065       ocoff(is)-=sum*ocoff(ik);
00066     }
00067     sum=norm2(ocoff(is));
00068     ocoff(is)=ocoff(is)/sqrt(sum);
00069   }
00070   return trans(ocoff.sub(skip,deg));
00071 }
00072 
00077 dmatrix orthpoly_constant_begin(int n,int deg,int nconst)
00078 {
00079   int j; int is; int ik;
00080   dmatrix ocoff(0,deg,1,n);
00081   double sum;
00082   ocoff(0)=sqrt(double(n));
00083   if (nconst>n-1)
00084   {
00085     cerr << "nconst too large in orthpoly_constant_begin"
00086          << endl;
00087   }
00088   if (deg>n-nconst)
00089   {
00090     cerr << "deg too large in orthpoly_constant_begin"
00091          << endl;
00092   }
00093   for (is=1; is<=deg; is++)
00094   {
00095     if (nconst>1)
00096     {
00097       for (j=1; j<=nconst; j++)
00098       {
00099         ocoff(is,j)=1.0;
00100       }
00101       for (j=nconst+1; j<=n; j++)
00102       {
00103         ocoff(is,j)=pow(double(j-nconst+1),is);
00104       }
00105     }
00106     else
00107     {
00108       for (j=1; j<=n; j++)
00109       {
00110         ocoff(is,j)=pow(double(j),is);
00111       }
00112     }
00113   }
00114   for (is=0; is<=deg; is++) /* L1000  */
00115   {
00116     for (ik=0; ik<=is-1; ik++) /* L2000  */
00117     {
00118       sum=ocoff(is)*ocoff(ik);
00119       ocoff(is)-=sum*ocoff(ik);
00120     }
00121     sum=norm2(ocoff(is));
00122     ocoff(is)=ocoff(is)/sqrt(sum);
00123   }
00124   int ps=0;
00125   if (ps)
00126   {
00127     dmatrix tmp(0,deg,0,deg);
00128     for (int i=0;i<=deg;i++)
00129     {
00130       for (int j=0;j<=deg;j++)
00131       {
00132         tmp(i,j)=ocoff(i)*ocoff(j);
00133       }
00134     }
00135     cout << tmp << endl;
00136   }
00137   return trans(ocoff);
00138 }
00139 
00144 dmatrix orthpoly_constant_begin_end(int n,int deg,int nconst_begin,
00145   int end_degree,int nconst_end)
00146 {
00147   int j; int is; int ik;
00148   dmatrix ocoff(0,deg,1,n);
00149   double sum;
00150   ocoff(0)=sqrt(double(n));
00151   if (nconst_begin>n-1)
00152   {
00153     cerr << "nconst_begin too large in orthpoly_constant_begin"
00154          << endl;
00155   }
00156   if (deg>n-nconst_begin)
00157   {
00158     cerr << "deg too large in orthpoly_constant_begin"
00159          << endl;
00160   }
00161   for (is=1; is<=deg; is++)
00162   {
00163     if (nconst_begin>1)
00164     {
00165       for (j=1; j<=nconst_begin; j++)
00166       {
00167         ocoff(is,j)=1.0;
00168       }
00169       for (j=nconst_begin+1; j<=n; j++)
00170       {
00171         int jj=j;
00172         if (j>n-nconst_end+1 && is>=end_degree)
00173         {
00174           jj=n-nconst_end+1;
00175         }
00176         ocoff(is,j)=pow(double(jj-nconst_begin+1)/n,is);
00177       }
00178     }
00179     else
00180     {
00181       for (j=1; j<=n; j++)
00182       {
00183         int jj=j;
00184         if (j>n-nconst_end+1 && is>=end_degree)
00185         {
00186           jj=n-nconst_end+1;
00187         }
00188         ocoff(is,j)=pow(double(jj)/n,is);
00189       }
00190     }
00191   }
00192   for (is=0; is<=deg; is++) /* L1000  */
00193   {
00194     for (ik=0; ik<=is-1; ik++) /* L2000  */
00195     {
00196       sum=ocoff(is)*ocoff(ik);
00197       ocoff(is)-=sum*ocoff(ik);
00198     }
00199     sum=norm2(ocoff(is));
00200     ocoff(is)=ocoff(is)/sqrt(sum);
00201   }
00202   int ps=0;
00203   if (ps)
00204   {
00205     dmatrix tmp(0,deg,0,deg);
00206     for (int i=0;i<=deg;i++)
00207     {
00208       for (int j=0;j<=deg;j++)
00209       {
00210         tmp(i,j)=ocoff(i)*ocoff(j);
00211       }
00212     }
00213     cout << tmp << endl;
00214   }
00215   return trans(ocoff);
00216 }
00217 
00222 dmatrix seldif_basis(int n)
00223 {
00224   int i; int j;
00225   dmatrix ocoff(1,n,1,n);
00226   dmatrix ocoff1(1,n,1,n);
00227   ocoff.initialize();
00228   ocoff1.initialize();
00229   for (i=1; i<=n; i++)
00230   {
00231     for (j=i; j<=n; j++)
00232     {
00233       ocoff(i,j)=1;
00234     }
00235   }
00236   ocoff1=trans(ocoff);
00237 
00238   for (i=1; i<=n; i++) /* L1000  */
00239   {
00240     for (j=1; j<=i-1; j++) /* L2000  */
00241     {
00242       ocoff(i)-=(ocoff(i)*ocoff(j))*ocoff(j);
00243     }
00244     ocoff(i)/=norm(ocoff(i));
00245   }
00246   ocoff=trans(ocoff);
00247 
00248   cout << setw(10) << setprecision(4) << ocoff1 << endl << endl;
00249 
00250   dmatrix tmp1=(inv(ocoff1)*ocoff);
00251   dvector a(1,n);
00252   dvector b(1,n);
00253 
00254   for (i=1; i<=n; i++) /* L1000  */
00255   {
00256     a(i)=tmp1(i,i);
00257   }
00258   b(1)=0.0;
00259   for (i=2; i<=n; i++) /* L1000  */
00260   {
00261     b(i)=tmp1(i-1,i);
00262   }
00263 
00264   cout << a << endl << endl;
00265   cout << b << endl << endl;
00266 
00267   cout << ocoff1*tmp1(1) << endl;
00268   cout << ocoff1*tmp1(2) << endl;
00269   cout << ocoff1 *tmp1(3) << endl;
00270   cout << (ocoff1*tmp1(1)) * (ocoff1 *tmp1(3)) << endl;
00271   cout << (ocoff1*tmp1(2)) * (ocoff1 *tmp1(3)) << endl;
00272 
00273   return ocoff1;
00274 }