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