Changeset bc1dd97 in tol

Timestamp:

Apr 24, 2012, 5:33:49 PM (15 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

1941-docker-image-to-test, devel, master

Children:

1466af1

Parents:

736be83

Message:

Refs #1511

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@4587 be4846ae-f051-483a-831f-b64e7d08527f

Location:

tol

Files:

• bmath/bstat/ar.cpp (modified) (1 diff)
• bmath/bstat/arma.cpp (modified) (5 diffs)
• bmath/bstat/tol_bar.h (modified) (1 diff)
• bmath/bstat/unitroot.cpp (modified) (2 diffs)
• bmath/mathobjects/tol_bpolyn.h (modified) (3 diffs)
• btol/bmodel/modcalc.cpp (modified) (4 diffs)
• btol/bmodel/moddiag.cpp (modified) (2 diffs)
• btol/real_type/datgrav.cpp (modified) (6 diffs)

tol/bmath/bstat/ar.cpp

@@ -298 +298 @@
     IdentifyAutoRegIntegrated(z,err,s,sw,w,dif,ar,period);
 }
+
+//--------------------------------------------------------------------
+BArray<BDat> NextSchur(const BArray<BDat>& pol)
+//--------------------------------------------------------------------
+{
+  int k, s = pol.Size();
+  BArray<BDat> b(s-1);
+  b.AllocBuffer(s-1);
+  double* B  = (double*)b.GetBuffer();
+  const double* A  = (const double*)pol.Buffer();
+  const double* A_ = (const double*)pol.Buffer()+s-1;
+  double a = *A; 
+  double a_ = *A_;
+  double c = a / a_;
+/* * /
+  Std(BText("\nNextSchur a0=")+a+" an="+a_+" c="+c);
+  if(BDat(c).IsUnknown())
+   printf("");
+/* */
+  A++;
+  A_--;
+  for(k=0; k<s-1; k++, A++, A_--, B++)
+  {
+    *B = (*A - c * *A_);
+  }
+  return(b);
+}
+
+//--------------------------------------------------------------------
+BDat SchurValue(const BArray<BDat>& a, BDat toleranceBase, BDat& tolerance, bool &isSchur)
+//--------------------------------------------------------------------
+{
+  int sld = sizeof(long double);
+  int s = a.Size();
+  BDat a0 = Abs(a[0]);
+  BDat an = Abs(a[s-1]);
+  BDat sv = (2*an)/(an+a0);
+  isSchur = an-a0 > tolerance;
+/* * /
+//Std(BText("\nSchurValue a=")); for(int k=0; k<s; k++) { Std(BText("|")+a[k]); }
+  Std(BText("\nSchurValue tolerance=")+tolerance+" a0=")+a0+" an="+an+" sv="+sv+" isSchur="+isSchur);
+/* */
+  if(!isSchur || (s==2))
+  {
+    return(sv);
+  } 
+  else 
+  {
+    tolerance += toleranceBase*(1+2*(s-1));
+    return(SchurValue(NextSchur(a),toleranceBase,tolerance,isSchur));
+  }
+}
+
+
+//--------------------------------------------------------------------
+BDat StationaryValue(
+  const BPolyn<BDat>& pol, 
+  BDat toleranceBase, 
+  BDat& tolerance, 
+  bool& isStationary)
+//--------------------------------------------------------------------
+{
+  
+  if(toleranceBase.IsUnknown() || toleranceBase<0) { toleranceBase = DEpsilon(); }
+  tolerance = toleranceBase;
+  int i,j, deg;
+  BDat coef;
+  int s = pol.Size();
+  int n = pol.Degree();
+  if(s==1)
+  {
+    if(n==0) { isStationary=true;  return(2); }
+    else     { isStationary=false; return(0); }
+  }
+  int gcd = pol.Period();
+  int m = n/gcd;
+  BArray<BDat> a(BDat(0),m+1);
+  for(i=0; i<s; i++)
+  {
+    j = s-1-i;
+    deg  = pol(i).Degree();
+    coef = pol(i).Coef();
+    a[m-deg/gcd] = coef;
+  }
+  BDat sv = SchurValue(a,toleranceBase,tolerance,isStationary);
+  if(sv.IsUnknown())
+  {
+    Warning("StationaryValue returns unknown value.");
+  }
+  return(sv);  
+}
+
+//--------------------------------------------------------------------
+bool IsStationary(const BPolyn<BDat>& pol, BDat& tolerance)
+//--------------------------------------------------------------------
+{
+  bool isStationary;
+  BDat toleranceBase = tolerance;
+  BDat sv = StationaryValue(pol,toleranceBase,tolerance,isStationary); 
+  return(isStationary);
+}
+
+//--------------------------------------------------------------------
+BDat StationaryValue(const BPolyn<BDat>& pol)
+//--------------------------------------------------------------------
+{
+  BDat tolerance, toleranceBase = BDat::Unknown();
+  bool isStationary;
+  return(StationaryValue(pol,toleranceBase,tolerance,isStationary));
+}
+
+//--------------------------------------------------------------------
+bool IsStationary(const BPolyn<BDat>& pol)
+//--------------------------------------------------------------------
+{
+  BDat tolerance = BDat::Unknown();
+  return(IsStationary(pol,tolerance));
+}

tol/bmath/bstat/arma.cpp

@@ -35 +35 @@
 #include <tol/tol_bvmat.h>
 #include <tol/tol_bratgra.h>
+#include <tol/tol_bar.h>
 
 //--------------------------------------------------------------------
@@ -330 +331 @@
   MatBackwardDifEq (R, ZEROp, ZB<<Z, ZEROq, A0A);
 
-  if(ar.IsStationary() && ma.IsStationary())
+  if(IsStationary(ar) && IsStationary(ma))
   {
     BMatrix<BDat> oldA0A;
@@ -700 +701 @@
   int i, j;
   BArray<BDat> gn;
-//BBool isStationary = ar.IsStationary() && ma.IsStationary();
+//BBool isStationary = IsStationary(ar) && IsStationary(ma);
   ARMAAutoCovarianzeVector(gn, ar, ma, n0, sigma);
   BSymMatrix<BDat> G(n0);
@@ -1283 +1284 @@
     BPolyn<BDat> ar_i = factor_(i).ar_; ar_i.Aggregate();
     BPolyn<BDat> ma_i = factor_(i).ma_; ma_i.Aggregate();
-    arOk = arOk && ar_i.IsStationary();
-    maOk = maOk && ma_i.IsStationary();
+    arOk = arOk && IsStationary(ar_i);
+    maOk = maOk && IsStationary(ma_i);
   }
   bool ok = arOk && maOk;
@@ -1297 +1298 @@
       ar_i.Aggregate();
       ma_i.Aggregate();
-      bool arOk = ar_i.IsStationary()!=0;
-      bool maOk = ma_i.IsStationary()!=0;
+      bool arOk = IsStationary(ar_i)!=0;
+      bool maOk = IsStationary(ma_i)!=0;
       Std(BText("\nStationary?[")<<(arOk?"YES":"NOT")<<"] AR("<<factor_(i).s_<<")="<<factor_(i).ar_.Name());
       Std(BText("\nStationary?[")<<(maOk?"YES":"NOT")<<"] MA("<<factor_(i).s_<<")="<<factor_(i).ma_.Name());

tol/bmath/bstat/tol_bar.h

@@ -67 +67 @@
                                BInt          period=1);
 
+BDat StationaryValue(
+  const BPolyn<BDat>& pol, 
+  BDat toleranceBase, 
+  BDat& tolerance, 
+  bool& isStationary);
+bool IsStationary(const BPolyn<BDat>& pol, BDat& tolerance);
+BDat StationaryValue(const BPolyn<BDat>& pol);
+bool IsStationary(const BPolyn<BDat>& pol);
+
+
 
 #endif  // TOL_BAR_H

tol/bmath/bstat/unitroot.cpp

@@ -104 +104 @@
     else
     {
-        ok = p.IsStationary();
+        ok = IsStationary(p);
     }
     if(changed && name.HasName())
@@ -180 +180 @@
         rAr[0].PutCoef(1);
         for(i=0; i<q-d; i++) { rAr[i+1].PutCoef(-(rNU(i,0)+NU(i,0))); }
-        if(rAr.IsStationary()) { prob+=1; }
+        if(IsStationary(rAr)) { prob+=1; }
     }
     detChoCOV = 1;

tol/bmath/mathobjects/tol_bpolyn.h

@@ -98 +98 @@
   void ChangeBF(const BPolyn<Any>& pol);
   void InverseRoots(const BPolyn<Any>& pol);
-  BPolyn<Any> NextSchur() const;
-  Any    Schur(BBool complete) const;
-  Any    StationaryValue(BBool complete) const;
-  BBool IsStationary() const;
 
   BBool IsKnown() const {
@@ -348 +344 @@
 }
 
-//--------------------------------------------------------------------
-template <class Any>
-BPolyn<Any> BPolyn<Any>::NextSchur() const
-//--------------------------------------------------------------------
-{
-//Std(BText("\nTRACE  NextSchur 1 this=")+Name());
-  BInt n = Degree();
-  BArray<Any> a(n+1); 
-  int k;
-  for(k=0; k<=n; k++) { a[k] = Coef(k); }
-  Any a0  = a[0];
-  Any an  = a[n];
-  Any c;
-//Std(BText("\na[0]: ")+a0.Name()+"; a[n]:"+an.Name());
-  BPolyn<Any> next;
-  next.ReallocBuffer(n);
-  for(k=0; k<n; k++)
-  {
-    c = an*a[k+1] - a0*a[n-k-1];
-    next(k).PutCoef(c);
-    next(k).PutDegree(k);
-  //Std(BText("\na[k+1]: ")+a[k+1].Name()+"; a[n-k-1]:"+a[n-k-1].Name()+"; next(i-1):"+c.Name());
-  }
-//Std(BText("\nTRACE NextSchur 2 next=")+next.Name());
-  next.Aggregate();
-//Std(BText("\nNextSchur 3 next=")+next.Name());
-//if(next.Degree()>0) { next /= next.Coef(next.Degree()); }
-//Std(BText("\nTRACE NextSchur 4 next=")+next.Name());
-  return(next);
-}
-
-
-//--------------------------------------------------------------------
-template <class Any>
-Any BPolyn<Any>::Schur(BBool complete) const
-
-/*! Returns a value great than 1 if there are no roots out of
- *  the unit circle
- */
-//--------------------------------------------------------------------
-{
-//Std(BText("\nTRACE  Schur : ")+Name()); 
-  int s = this->Size();
-  int n = this->Degree();
-  if(s==1) { return(2); }
-  Any a0 = Abs(this->Coef(0));
-  Any an = Abs(this->Coef(n));
-  Any v = 1.0+(an-a0)/(an+a0);
-  if((!complete && (v<=1)) || (s==2))
-  {
-    return(v);
-  } 
-  else 
-  {
-    return(NextSchur().Schur(complete));
-  }
-}
 
 //--------------------------------------------------------------------
@@ -423 +362 @@
 }
 
-//--------------------------------------------------------------------
-template <class Any>
-Any BPolyn<Any>::StationaryValue(BBool complete) const
-
-/*! Returns a value great than 1 if there are no roots inside
- *  the unit circle
- */
-//--------------------------------------------------------------------
-{
-//Std(BText("\nTRACE  StationaryValue : ")+Name());
-  int s = this->Size();
-  int n = this->Degree();
-  if(s==1)
-  {
-    if(n==0) { return(2); }
-    else     { return(0); }
-  }
-  int i, j;
-  int deg;
-  Any coef;
-  BPolyn sym;
-  int gcd = Period();
-//Std(BText("\nTRACE  StationaryValue  gcd=")+gcd);
-  int m = n/gcd;
-  sym.ReallocBuffer(s);
-  for(i=0; i<s; i++)
-  {
-    j = s-1-i;
-    deg  = (*this)(i).Degree();
-    coef = (*this)(i).Coef();
-    sym(j).PutDegree(m-deg/gcd);
-    sym(j).PutCoef(coef);
-  }
-  return(sym.Schur(complete));
-}
-
-
-//--------------------------------------------------------------------
-template <class Any>
-BBool BPolyn<Any>::IsStationary() const
-
-/*! Returns true if all roots in this polyn are outside of the unit
- *  circle
- */
-//--------------------------------------------------------------------
-{
-  if(this->Size  ()==0) { return(BFALSE); }
-  if(this->Degree()==0) { return(BTRUE); }
-  else { return(this->StationaryValue(BFALSE)>1); }
-}
-
 
 //--------------------------------------------------------------------

tol/btol/bmodel/modcalc.cpp

@@ -33 +33 @@
 #include <tol/tol_bsparse.h>
 #include <tol/tol_bsvdsep.h>
+#include <tol/tol_bar.h>
 
 #include <gsl/gsl_math.h>
@@ -1003 +1004 @@
   for(i=0; i<arma_.factor_.Size(); i++)
   {
-    bool arOk = arma_.factor_(i).ar_.IsStationary();
-    bool maOk = arma_.factor_(i).ma_.IsStationary();
+    bool arOk = IsStationary(arma_.factor_(i).ar_);
+    bool maOk = IsStationary(arma_.factor_(i).ma_);
     Std(BText("\nStationary?[")<<(arOk?"YES":"NOT")<<"] AR("<<arma_.factor_(i).s_<<")="<<arma_.factor_(i).ar_.Name());
     Std(BText("\nStationary?[")<<(maOk?"YES":"NOT")<<"] MA("<<arma_.factor_(i).s_<<")="<<arma_.factor_(i).ma_.Name());
@@ -1104 +1105 @@
   for(i=0; i<arma_.factor_.Size(); i++)
   {
-    bool arOk = arma_.factor_(i).ar_.IsStationary();
-    bool maOk = arma_.factor_(i).ma_.IsStationary();
+    bool arOk = IsStationary(arma_.factor_(i).ar_);
+    bool maOk = IsStationary(arma_.factor_(i).ma_);
     Std(BText("\nStationary?[")<<(arOk?"YES":"NOT")<<"] AR("<<arma_.factor_(i).s_<<")="<<arma_.factor_(i).ar_.Name());
     Std(BText("\nStationary?[")<<(maOk?"YES":"NOT")<<"] MA("<<arma_.factor_(i).s_<<")="<<arma_.factor_(i).ma_.Name());
@@ -1297 +1298 @@
 //--------------------------------------------------------------------
 {
-  BDat y = p_.StationaryValue(BTRUE);
+  BDat y = StationaryValue(p_);
 //Std(BText("\nS(")+p_.Name()+")="+(y-1));
   return(y-1);

tol/btol/bmodel/moddiag.cpp

@@ -33 +33 @@
 #include <tol/tol_bmultvar.h>
 #include <tol/tol_bsetgra.h>
+#include <tol/tol_bar.h>
 
 //--------------------------------------------------------------------
@@ -468 +469 @@
 //        Std(BText("\nFactor[")+j+"]"+factors[j].Name() +
 //                " Stationary = " + factors[j].Stationary());
-        if(!(factors[j].IsStationary())) { unarys[j]+=1; }
+        if(!(IsStationary(factors[j]))) { unarys[j]+=1; }
       }
     }

tol/btol/real_type/datgrav.cpp

@@ -40 +40 @@
 #include <tol/tol_bdtegra.h>
 #include <tol/tol_bnameblock.h>
+#include <tol/tol_bar.h>
 #ifdef __USE_TC__
 #  include <tol/tol_bctmigra.h>
@@ -1311 +1312 @@
 //--------------------------------------------------------------------
   DeclareContensClass(BDat, BDatTemporary, BDatIsStationary);
-  DefExtOpr(1, BDatIsStationary, "IsStationary", 1, 1, "Polyn",
-  "(Polyn pol)",
+  DefExtOpr(1, BDatIsStationary, "IsStationary", 1, 2, "Polyn Real",
+  "(Polyn pol [, Real tolerance=MachineEpsilon])",
   I2("Returns true if all the roots of the polynomial are would be of the "
      "unit circle.",
@@ -1321 +1322 @@
 //--------------------------------------------------------------------
 {
-  BPol p = Pol(Arg(1));
-  p.Aggregate();
-  contens_ = p.IsStationary();
+  BPol& p = Pol(Arg(1));
+  BDat t;
+  if(Arg(2)) { t = Dat(Arg(2)); }
+  contens_ = IsStationary(p,t);
 }
 
@@ -1339 +1341 @@
 //--------------------------------------------------------------------
 {
-  contens_ = Pol(Arg(1)).Period();
+  BPol& p = Pol(Arg(1));
+  contens_ = p.Period();
 }
 
 //--------------------------------------------------------------------
   DeclareContensClass(BDat, BDatTemporary, BDatStationaryValue);
-  DefExtOpr(1, BDatStationaryValue, "StationaryValue", 1, 1, "Polyn",
-  "(Polyn pol)",
+  DefExtOpr(1, BDatStationaryValue, "StationaryValue", 1, 2, "Polyn Real",
+  "(Polyn pol [, Real tolerance=MachineEpsilon])",
   I2("Returns a measure of the distence to the unit circle of the "
      "polynomial roots .",
@@ -1354 +1357 @@
 //--------------------------------------------------------------------
 {
-  contens_ = Pol(Arg(1)).StationaryValue(BFALSE);
+  BPol& p = Pol(Arg(1));
+  BDat toleranceBase;
+  BDat tolerance;
+  if(Arg(2)) { toleranceBase = Dat(Arg(2)); }
+  bool isStationary;
+  contens_ = StationaryValue(p,toleranceBase,tolerance,isStationary);
 }
 
@@ -1413 +1421 @@
     }
     pk.Aggregate();
-    BInt ok = pk.IsStationary();
+    BInt ok = IsStationary(pk);
 //  Std(BText("\n")+(ok?" OK   ":" FAIL ")+pk.Name());
     sumOk += ok;