Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members  

SparseVec.cc

Go to the documentation of this file.
00001 /*
00002     File:           SparseVec.cc
00003 
00004     Function:       Implements SparseVec.h
00005 
00006     Author(s):      Andrew Willmott
00007 
00008     Copyright:      (c) 1995-2000, Andrew Willmott
00009 
00010     Notes:          
00011 
00012 */
00013 
00014 
00015 #include "vl/SparseVec.h"
00016 #include <stdarg.h>
00017 #include <iomanip.h>
00018 #include "vl/SparseVec.h"
00019 #include "vl/SubSVec.h"
00020 
00021 
00022 TSparseVec::TSparseVec() : TSVList(), elts(0)
00023 {
00024     End();
00025 }
00026 
00027 TSparseVec::TSparseVec(Int n) : TSVList(), elts(n)
00028 {
00029     Assert(n > 0,"(SparseVec) illegal vector size");
00030 }
00031 
00032 TSparseVec::TSparseVec(const TSparseVec &v) : TSVList(v), elts(v.elts)
00033 {
00034 }
00035 
00036 TSparseVec::TSparseVec(const TSubSVec &v) : TSVList()
00037 {
00038     SELF = v;
00039 }
00040 
00041 TSparseVec::TSparseVec(const TVec &v) : TSVList()
00042 {
00043     SELF = v;
00044 }
00045 
00046 TSparseVec::TSparseVec(Int n, Int indices[], TVReal elts[]) : TSVList(), elts(n)
00047 {
00048     Assert(n > 0,"(SparseVec) illegal vector size");
00049     Int i = 0;
00050     
00051     while (1)
00052     {
00053         if (indices[i] < 0)
00054             break;
00055         AddElt(indices[i], elts[i]);
00056         i++;
00057     }
00058     
00059     End();
00060 }
00061 
00062 TSparseVec::TSparseVec(Int n, Int idx0, double elt0, ...) : TSVList(), elts(n)
00063 {
00064     Assert(n > 0,"(SparseVec) illegal vector size");
00065     va_list     ap;
00066     Int         idx;
00067     TVReal      elt;
00068     
00069     va_start(ap, elt0);
00070     Assert(idx0 >= 0 && idx0 < elts, "(SparseVec) illegal first index");
00071     SetReal(elt, elt0);
00072     AddElt(idx0, elt);
00073         
00074     while (1)
00075     {
00076         idx = va_arg(ap, Int);
00077         if (idx < 0)
00078             break;
00079         Assert(idx < elts, "(SparseVec) illegal index");
00080         SetReal(elt, va_arg(ap, double));
00081         AddElt(idx, elt);
00082     }
00083     
00084     End();
00085     va_end(ap);
00086 }
00087 
00088 TSparseVec::TSparseVec(Int n, ZeroOrOne k) : elts(n)
00089 {
00090     Assert(n > 0,"(SparseVec) illegal vector size");
00091     MakeBlock(k);
00092 }
00093 
00094 TSparseVec::TSparseVec(Int n, Axis a) : elts(n)
00095 {
00096     Assert(n > 0,"(SparseVec) illegal vector size");
00097     MakeUnit(a);
00098 }
00099 
00100 
00101 TSparseVec::~TSparseVec()
00102 {
00103 }
00104 
00105 TSparseVec &TSparseVec::operator = (const TSparseVec &v)
00106 {
00107     Int i;
00108     
00109     SetNumElts(v.Elts());
00110     Begin();
00111     
00112     for (i = 0; i < v.NumItems() - 1; i++)
00113         Append(v[i]);
00114     
00115     End();
00116     
00117     return(SELF);
00118 }
00119 
00120 TSparseVec &TSparseVec::operator = (const TVec &v)
00121 {
00122     Int i;
00123 
00124     SetNumElts(v.Elts());
00125     Begin();
00126     
00127     for (i = 0; i < v.Elts(); i++)  
00128         AddElt(i, v[i]);
00129     
00130     End();
00131     
00132     return(SELF);
00133 }
00134 
00135 TSparseVec &TSparseVec::operator = (const TSubSVec &v)
00136 {
00137     v.Store(SELF);
00138     return(SELF);
00139 }
00140 
00141 Void TSparseVec::SetSize(Int n)
00142 {
00143     SetNumElts(n);
00144     Begin();    
00145     End();
00146 }
00147 
00148 Void TSparseVec::MakeZero()
00149 {
00150     Begin();
00151     End();
00152 }
00153 
00154 Void TSparseVec::MakeUnit(Int i, TVReal k)
00155 {
00156     Begin();
00157     AddElt(i, k);
00158     End();
00159 }
00160 
00161 Void TSparseVec::MakeBlock(TVReal k)
00162 {
00163     if (len(k) == 0.0)
00164         MakeZero();
00165     else
00166         _Error("(SparseVec::MakeBlock) Calling this method is a really dumb idea.");
00167     // What's the point of setting every member of a sparse vector to
00168     // something non-zero? 
00169 }
00170 
00171 
00172 Bool TSparseVec::compactPrint = false;
00173 TVReal TSparseVec::fuzz = 1e-10;
00174 
00175 Void TSparseVec::SetCompactPrint(Bool on)
00176 { compactPrint = on; }
00177 
00178 Void TSparseVec::SetFuzz(TVReal theFuzz)
00179 { fuzz = theFuzz; }
00180 
00181 
00182 // --- Vector operations ------------------------------------------------------
00183 
00184 
00185 #include <ctype.h>
00186 
00187 /*  [Note]
00188 
00189     When operating on sparse vectors in place, it is faster to create a 
00190     new vector on the fly, and swap it in at the end, than to perform
00191     inserts/deletes on the original. (O(n) vs O(n^2).)
00192     
00193     Supporting Analysis:
00194     
00195     There will be O(n) new/deleted elts. 
00196     Copying will involve O(n) copies in creating the new array, plus 
00197     a temporary memory overhead the size of the old vector. Operating
00198     in place will require O(n) insert/deletes, each of which require
00199     O(n) copies, and thus requires O(n^2) copies overall.
00200 */
00201 
00202 TSparseVec &operator += (TSparseVec &a, const TSparseVec &b)
00203 {
00204     TSparseVec  t = a + b;  
00205     a.Replace(t);
00206     return(a);
00207 }
00208 
00209 TSparseVec &operator -= (TSparseVec &a, const TSparseVec &b)
00210 {
00211     TSparseVec  t = a - b;
00212     a.Replace(t);
00213     return(a);
00214 }
00215 
00216 TSparseVec &operator *= (TSparseVec &a, const TSparseVec &b)
00217 {
00218     TSparseVec  t = a * b;
00219     a.Replace(t);
00220     return(a);
00221 }
00222 
00223 TSparseVec &operator *= (TSparseVec &v, TVReal s)
00224 {
00225     TSparseVec  t = v * s;
00226     v.Replace(t);
00227     return(v);
00228 }
00229 
00230 TSparseVec &operator /= (TSparseVec &a, const TSparseVec &b)
00231 {
00232     TSparseVec  t = a / b;
00233     a.Replace(t);
00234     return(a);
00235 }
00236 
00237 TSparseVec &operator /= (TSparseVec &v, TVReal s)
00238 {
00239     TSparseVec t = v / s;
00240     v.Replace(t);
00241     return(v);
00242 }
00243 
00244 
00245 // --- Vec Comparison Operators -----------------------------------------------
00246 
00247 
00248 Bool operator == (const TSparseVec &a, const TSparseVec &b)
00249 {
00250     Assert(a.Elts() == b.Elts(), "(SparseVec::==) Vec sizes don't match");
00251 
00252     if (a.NumItems() != b.NumItems())
00253         return(false);
00254     else
00255     {
00256         Int i;
00257         
00258         for (i = 0; i < a.NumItems(); i++)
00259             if (a[i] != b[i])
00260                 return(false);
00261         
00262         return(true);
00263     }
00264 }
00265 
00266 Bool operator != (const TSparseVec &a, const TSparseVec &b)
00267 {
00268     return(!(a == b));
00269 }
00270 
00271 
00272 // --- SparseVec Methods ------------------------------------------------------
00273 
00274 
00275 TSparseVec TSparseVec::Overlay(const TSparseVec &b) const
00276 {
00277     TSparseVec  result(Elts());
00278     Int         i, j;
00279 
00280     for (i = 0, j = 0; ; )
00281         if (item[i].index == b[j].index)
00282         {
00283             if (item[i].index == VL_SV_MAX_INDEX)
00284                 break;
00285             
00286             result.AddNZElt(b[j]);
00287             i++;
00288             j++;
00289         }
00290         else if (item[i].index < b[j].index)
00291         {   result.AddNZElt(item[i]); i++;  }
00292         else
00293         {   result.AddNZElt(b[j]); j++; }
00294             
00295     result.End();   
00296     return(result);
00297 }
00298 
00299 Void TSparseVec::SetElts(Int idx0, double elt0, ...)
00300 {
00301     va_list     ap;
00302     Int         idx;
00303     TVReal      elt;
00304     
00305     va_start(ap, elt0);
00306     Begin();
00307     Assert(idx0 >= 0 && idx0 < elts, "(SparseVec::SetElts) illegal first index");
00308     SetReal(elt, elt0);
00309     AddElt(idx0, elt);
00310         
00311     while (1)
00312     {
00313         idx = va_arg(ap, Int);
00314         if (idx < 0)
00315             break;
00316         Assert(idx < elts, "(SparseVec::SetElts) illegal index");
00317         SetReal(elt, va_arg(ap, double));
00318         AddElt(idx, elt);
00319     }
00320     
00321     End();
00322     va_end(ap);
00323 }
00324 
00325 Void TSparseVec::Set(Int index, TVReal elt)
00326 {
00327     TSVIter     j(SELF);
00328     
00329     if (len(elt) <= fuzz)
00330         return;
00331                     
00332     if (!j.IncTo(index))
00333     {
00334         Insert(j.PairIdx());
00335         Item(j.PairIdx()).index = index;
00336     }
00337 
00338     Item(j.PairIdx()).elt = elt;
00339 }
00340 
00341 TVReal TSparseVec::Get(Int index) const
00342 {
00343     TSVIter     j(SELF);
00344     
00345     if (j.IncTo(index))
00346         return(j.Data());
00347     else
00348         return(vl_zero);
00349 }
00350 
00351 // --- Vec Arithmetic Operators -----------------------------------------------
00352 
00353 TSparseVec operator + (const TSparseVec &a, const TSparseVec &b)
00354 {
00355     Assert(a.Elts() == b.Elts(), "(SparseVec::+) Vec sizes don't match");
00356 
00357     TSparseVec  result(a.Elts());
00358     Int         i, j;
00359 
00360     // Step through a and b in parallel
00361 
00362     for (i = 0, j = 0; ; )
00363         if (a[i].index == b[j].index)
00364         {
00365             // We have two elements at the same index. 
00366             // Are we at the end of both arrays?
00367             if (a[i].index == VL_SV_MAX_INDEX)
00368                 break;
00369             
00370             // If not, add the result
00371             
00372             result.AddElt(a[i].index, a[i].elt + b[j].elt); // +
00373             i++;
00374             j++;
00375         }
00376         else if (a[i].index < b[j].index)
00377         // result[x] = a[i] + 0     
00378         {   result.AddNZElt(a[i]); i++; }
00379         else
00380         // result[x] = b[j] + 0     
00381         {   result.AddNZElt(b[j]); j++; }
00382             
00383     result.End();   
00384     return(result);
00385 }
00386 
00387 TSparseVec operator - (const TSparseVec &a, const TSparseVec &b) 
00388 {
00389     Assert(a.Elts() == b.Elts(), "(SparseVec::-) Vec sizes don't match");
00390 
00391     TSparseVec  result(a.Elts());
00392     Int         i, j;
00393 
00394     for (i = 0, j = 0; ; )
00395         if (a[i].index == b[j].index)
00396         {
00397             if (a[i].index == VL_SV_MAX_INDEX)
00398                 break;
00399             
00400             result.AddElt(a[i].index, a[i].elt - b[j].elt); // -
00401             i++;
00402             j++;
00403         }
00404         else if (a[i].index < b[j].index)
00405         {   result.AddNZElt(a[i]); i++; }
00406         else
00407         {   result.AddNZElt(b[j].index, -b[j].elt); j++;    }
00408     
00409     result.End();   
00410     return(result);
00411 }
00412 
00413 TSparseVec operator - (const TSparseVec &v)
00414 {
00415     TSparseVec  result(v.Elts());
00416     Int         i;
00417     
00418     for (i = 0; i < v.NumItems() - 1; i++)
00419         result.AddNZElt(v[i].index, -v[i].elt);
00420     
00421     result.End();   
00422     return(result);
00423 }
00424 
00425 TSparseVec operator * (const TSparseVec &a, const TSparseVec &b)            
00426 {
00427     Assert(a.Elts() == b.Elts(), "(SparseVec::*) Vec sizes don't match");
00428 
00429     TSparseVec  result(a.Elts());
00430     TSVIter     j(a);
00431     Int         i;
00432     
00433     for (i = 0; i < b.NumItems() - 1; i++)
00434         if (j.IncTo(b[i].index))
00435             result.AddElt(b[i].index, j.Data() * b[i].elt);
00436     
00437     result.End();   
00438     return(result);
00439 }
00440 
00441 TSparseVec operator * (const TSparseVec &v, TVReal s) 
00442 {
00443     TSparseVec  result(v.Elts());
00444     Int         i;
00445         
00446     for (i = 0; i < v.NumItems() - 1; i++)
00447         result.AddElt(v[i].index, s * v[i].elt);
00448     
00449     result.End();   
00450     return(result);
00451 }
00452 
00453 TSparseVec operator / (const TSparseVec &a, const TSparseVec &b)            
00454 {
00455     Assert(a.Elts() == b.Elts(), "(SparseVec::/) Vec sizes don't match");
00456     
00457     TSparseVec  result(a.Elts());
00458     TSVIter     j(a);
00459     Int         i;
00460     
00461     for (i = 0; i < b.NumItems() - 1; i++)
00462         if (j.IncTo(b[i].index))
00463             result.AddElt(b[i].index, j.Data() / b[i].elt);
00464     
00465     result.End();   
00466     return(result);
00467 }
00468 
00469 TSparseVec operator / (const TSparseVec &v, TVReal s) 
00470 {
00471     TSparseVec  result(v.Elts());
00472     Int         i;
00473     TVReal      t = vl_1 / s;
00474             
00475     for (i = 0; i < v.NumItems() - 1; i++)
00476         result.AddElt(v[i].index, v[i].elt * t);
00477     
00478     result.End();   
00479     return(result);
00480 }
00481 
00482 TVReal dot(const TSparseVec &a, const TSparseVec &b) 
00483 {
00484     Assert(a.Elts() == b.Elts(), "(SparseVec::dot) Vec sizes don't match");
00485 
00486     TMReal      sum = vl_zero;
00487     TSVIter     j(a);
00488     Int         i;
00489     
00490     for (i = 0; i < b.NumItems() - 1; i++)
00491         if (j.IncTo(b[i].index))
00492             sum += j.Data() * b[i].elt;
00493     
00494     return(sum);
00495 }
00496 
00497 TVReal dot(const TSparseVec &a, const TVec &b)
00498 {
00499     Assert(a.Elts() == b.Elts(), "(SparseVec::dot) Vec sizes don't match");
00500 
00501     TMReal      sum = vl_zero;
00502     Int         i;
00503     
00504     for (i = 0; i < a.NumItems() - 1; i++)
00505         sum += a[i].elt * b[a[i].index];
00506     
00507     return(sum);
00508 }
00509 
00510 TSparseVec operator * (TVReal s, const TSparseVec &v)
00511 {
00512     return(v * s);
00513 }
00514 
00515 TVReal sqrlen(const TSparseVec &v)
00516 {
00517     TVReal      sum = vl_zero;
00518     Int         i;
00519     
00520     for (i = 0; i < v.NumItems() - 1; i++)
00521         sum += sqrlen(v[i].elt);
00522     
00523     return(sum);
00524 }
00525 
00526 TVReal len(const TSparseVec &v)
00527 {
00528     return(sqrt(sqrlen(v)));
00529 }
00530 
00531 TSparseVec norm(const TSparseVec &v)    
00532 {
00533     return(v / len(v));
00534 }
00535 
00536 
00537 // --- Vec Input & Output -----------------------------------------------------
00538 
00539 
00540 ostream &operator << (ostream &s, const TSparsePair &sp)
00541 {
00542     s << sp.index << ':' << sp.elt;
00543 
00544     return(s);
00545 }
00546 
00547 ostream &operator << (ostream &s, const TSparseVec &v)
00548 {
00549     if (TSparseVec::IsCompact())
00550     {
00551         Int i;
00552 
00553         s << '[';
00554 
00555         for (i = 0; i < v.NumItems() - 2; i++)
00556             s << v[i] << ' ';
00557 
00558         s << v[i] << ']';
00559     }
00560     else
00561     {
00562         Int i, j;
00563 
00564         s << '[';
00565         
00566         for (i = 0, j = 0; i < v.Elts() - 1; i++)
00567             if (i < v[j].index)
00568                 s << "0 ";
00569             else
00570                 s << v[j++].elt << ' ';
00571 
00572         if (i < v[j].index)
00573             s << "0]";
00574         else
00575             s << v[j].elt << ']';
00576     }
00577 
00578     return(s);
00579 }
00580 
00581 istream &operator >> (istream &s, TSparseVec &v)
00582 {
00583     Char    c;
00584     Int     i = 0;
00585     TVReal  elt;
00586 
00587     //  Expected format: [a b c d ...]
00588     
00589     while (isspace(s.peek()))           //  chomp white space
00590         s.get(c);
00591         
00592     if (s.peek() == '[')                        
00593     {
00594         v.Begin();
00595 
00596         s.get(c);
00597         
00598         while (isspace(s.peek()))       //  chomp white space
00599             s.get(c);
00600         
00601         while (s.peek() != ']')
00602         {           
00603             s >> elt;
00604             
00605             if (!s)
00606             {
00607                 _Warning("Couldn't read array component");
00608                 return(s);
00609             }
00610             else
00611             {
00612                 v.AddElt(i, elt);
00613                 i++;
00614             }
00615                 
00616             while (isspace(s.peek()))   //  chomp white space
00617                 s.get(c);
00618         }
00619         s.get(c);
00620     }
00621     else
00622     {
00623         s.clear(ios::failbit);
00624         _Warning("Error: Expected '[' while reading array");
00625         return(s);
00626     }
00627     
00628     v.End();
00629     v.SetNumElts(i);
00630     
00631     return(s);
00632 }
00633 
00634 
00635 // --- SparseVec iterator -----------------------------------------------------
00636 
00637 
00638 Void TSVIter::OrdFindLinear(Int i)
00639 {
00640     // Linear search for the right pair
00641         
00642     while (!AtEnd() && (i > Index()))
00643         pairIdx++;
00644 }
00645 
00646 #define SV_MIXED_SEARCH
00647 
00648 #ifndef __SV_CONST__
00649 #define __SV_CONST__
00650 static const Int kLinearSearchRange = 10;   //  Linear search over intervals smaller than this...
00651 static const Int kSparsenessEst = 5;        // estimated elts per non-zero elt.
00652 static const Int kLSRSparse = kSparsenessEst * kLinearSearchRange;
00653 #endif
00654 
00655 Void TSVIter::OrdFindBinary(Int i)
00656 //  Move index to point to the pair that corresponds to i, or if no such pair exists,
00657 //  the pair with the next index after i that does exist.
00658 {
00659 #ifdef SV_MIXED_SEARCH  
00660     // Mixture of linear & binary, parameterised by kLinearSearchRange.
00661     // If the item we're looking for is farther away from the current 
00662     // pair than kLSRSparse, we binary search.
00663 
00664     if ((i - Index()) > kLSRSparse)
00665     {
00666 #endif
00667 
00668     // --- Binary search on the pairs list-------------------------------------
00669     
00670     // A really nice thing to do here would be to back out
00671     // hierarchically from the current index instead of just
00672     // giving up on it. Similar to storing the current octree
00673     // node in RT, and doing the same.
00674 
00675         Int j = 0, k = list->NumItems() - 1, m;
00676         
00677         //  Test for trivial cases: i lies before or after the pairs list
00678         
00679         if (k < 0 || i <= list->Item(j).index)
00680         {   pairIdx = 0; return; }
00681         if (list->Item(k).index < i)
00682         {   pairIdx = k + 1; return; }
00683         if (list->Item(k).index == i)
00684         {   pairIdx = k; return; }
00685         
00686         while(k > j + 1 + kLinearSearchRange)
00687         {
00688             // precondition: j.index < i < k.index
00689             Assert(list->Item(j).index < i && i < list->Item(k).index, "Precondition failed.");
00690             
00691             // Naive midpoint picking
00692             m = (j + k) / 2;    // m ~ [j+1..k-1]
00693 
00694             // Linear midpoint picking
00695             // m = j + 1 +  ((k - j - 1) * (i - list->Item(j))) / list->Item(k) - list->Item(j));
00696                         
00697             if (i > list->Item(m).index)
00698                 j = m;
00699             else
00700             {
00701                 k = m;  
00702 
00703                 if (i >= list->Item(k).index)
00704                 {   
00705                     pairIdx = k;
00706                     return; 
00707                 }                   
00708             }
00709         }
00710 
00711         pairIdx = j + 1;
00712 
00713 #ifdef SV_MIXED_SEARCH
00714     }
00715 
00716     while (i > Index())     // Linear search, assuming sentinel
00717         pairIdx++;
00718 #endif
00719 }   
00720 

Generated at Sat Aug 5 00:16:48 2000 for Class Library by doxygen 1.1.0 written by Dimitri van Heesch, © 1997-2000