ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
orthpoly.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: orthpoly.cpp 2585 2014-11-07 19:12:51Z 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   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   int ps=0;
00122   if (ps)
00123   {
00124     dmatrix tmp(0,deg,0,deg);
00125     for (int i=0;i<=deg;i++)
00126     {
00127       for (int j=0;j<=deg;j++)
00128       {
00129         tmp(i,j)=ocoff(i)*ocoff(j);
00130       }
00131     }
00132     cout << tmp << endl;
00133   }
00134   return trans(ocoff);
00135 }
00136 
00141 dmatrix orthpoly_constant_begin_end(int n,int deg,int nconst_begin,
00142   int end_degree,int nconst_end)
00143 {
00144   dmatrix ocoff(0,deg,1,n);
00145   double sum;
00146   ocoff(0)=sqrt(double(n));
00147   if (nconst_begin>n-1)
00148   {
00149     cerr << "nconst_begin too large in orthpoly_constant_begin"
00150          << endl;
00151   }
00152   if (deg>n-nconst_begin)
00153   {
00154     cerr << "deg too large in orthpoly_constant_begin"
00155          << endl;
00156   }
00157   for (int is=1; is<=deg; is++)
00158   {
00159     if (nconst_begin>1)
00160     {
00161       for (int j=1; j<=nconst_begin; j++)
00162       {
00163         ocoff(is,j)=1.0;
00164       }
00165       for (int j=nconst_begin+1; j<=n; j++)
00166       {
00167         int jj=j;
00168         if (j>n-nconst_end+1 && is>=end_degree)
00169         {
00170           jj=n-nconst_end+1;
00171         }
00172         ocoff(is,j)=pow(double(jj-nconst_begin+1)/n,is);
00173       }
00174     }
00175     else
00176     {
00177       for (int j=1; j<=n; j++)
00178       {
00179         int jj=j;
00180         if (j>n-nconst_end+1 && is>=end_degree)
00181         {
00182           jj=n-nconst_end+1;
00183         }
00184         ocoff(is,j)=pow(double(jj)/n,is);
00185       }
00186     }
00187   }
00188   for (int is=0; is<=deg; is++) /* L1000  */
00189   {
00190     for (int ik=0; ik<=is-1; ik++) /* L2000  */
00191     {
00192       sum=ocoff(is)*ocoff(ik);
00193       ocoff(is)-=sum*ocoff(ik);
00194     }
00195     sum=norm2(ocoff(is));
00196     ocoff(is)=ocoff(is)/sqrt(sum);
00197   }
00198   int ps=0;
00199   if (ps)
00200   {
00201     dmatrix tmp(0,deg,0,deg);
00202     for (int i=0;i<=deg;i++)
00203     {
00204       for (int j=0;j<=deg;j++)
00205       {
00206         tmp(i,j)=ocoff(i)*ocoff(j);
00207       }
00208     }
00209     cout << tmp << endl;
00210   }
00211   return trans(ocoff);
00212 }
00213 
00218 dmatrix seldif_basis(int n)
00219 {
00220   dmatrix ocoff(1,n,1,n);
00221   dmatrix ocoff1(1,n,1,n);
00222   ocoff.initialize();
00223   ocoff1.initialize();
00224   for (int i=1; i<=n; i++)
00225   {
00226     for (int j=i; j<=n; j++)
00227     {
00228       ocoff(i,j)=1;
00229     }
00230   }
00231   ocoff1=trans(ocoff);
00232 
00233   for (int i=1; i<=n; i++) /* L1000  */
00234   {
00235     for (int j=1; j<=i-1; j++) /* L2000  */
00236     {
00237       ocoff(i)-=(ocoff(i)*ocoff(j))*ocoff(j);
00238     }
00239     ocoff(i)/=norm(ocoff(i));
00240   }
00241   ocoff=trans(ocoff);
00242 
00243   cout << setw(10) << setprecision(4) << ocoff1 << endl << endl;
00244 
00245   dmatrix tmp1=(inv(ocoff1)*ocoff);
00246   dvector a(1,n);
00247   dvector b(1,n);
00248 
00249   for (int i=1; i<=n; i++) /* L1000  */
00250   {
00251     a(i)=tmp1(i,i);
00252   }
00253   b(1)=0.0;
00254   for (int i=2; i<=n; i++) /* L1000  */
00255   {
00256     b(i)=tmp1(i-1,i);
00257   }
00258 
00259   cout << a << endl << endl;
00260   cout << b << endl << endl;
00261 
00262   cout << ocoff1*tmp1(1) << endl;
00263   cout << ocoff1*tmp1(2) << endl;
00264   cout << ocoff1 *tmp1(3) << endl;
00265   cout << (ocoff1*tmp1(1)) * (ocoff1 *tmp1(3)) << endl;
00266   cout << (ocoff1*tmp1(2)) * (ocoff1 *tmp1(3)) << endl;
00267 
00268   return ocoff1;
00269 }