1 /* 2 3 HyPhy - Hypothesis Testing Using Phylogenies. 4 5 Copyright (C) 1997-now 6 Core Developers: 7 Sergei L Kosakovsky Pond (spond@ucsd.edu) 8 Art FY Poon (apoon42@uwo.ca) 9 Steven Weaver (sweaver@ucsd.edu) 10 11 Module Developers: 12 Lance Hepler (nlhepler@gmail.com) 13 Martin Smith (martin.audacis@gmail.com) 14 15 Significant contributions from: 16 Spencer V Muse (muse@stat.ncsu.edu) 17 Simon DW Frost (sdf22@cam.ac.uk) 18 19 Permission is hereby granted, free of charge, to any person obtaining a 20 copy of this software and associated documentation files (the 21 "Software"), to deal in the Software without restriction, including 22 without limitation the rights to use, copy, modify, merge, publish, 23 distribute, sublicense, and/or sell copies of the Software, and to 24 permit persons to whom the Software is furnished to do so, subject to 25 the following conditions: 26 27 The above copyright notice and this permission notice shall be included 28 in all copies or substantial portions of the Software. 29 30 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 31 OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 32 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 33 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY 34 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, 35 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE 36 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 37 38 */ 39 40 #ifndef __MATRIX__ 41 #define __MATRIX__ 42 43 #include "hy_strings.h" 44 #include "avllistx.h" 45 #include "variablecontainer.h" 46 #include "trie.h" 47 48 #define _POLYNOMIAL_TYPE 0 49 #define _NUMERICAL_TYPE 1 50 #define _FORMULA_TYPE 2 51 #define _SIMPLE_FORMULA_TYPE 3 52 53 54 55 //_____________________________________________________________________________________________ 56 57 #define _HY_MATRIX_RANDOM_DIRICHLET 01L 58 #define _HY_MATRIX_RANDOM_GAUSSIAN 02L 59 #define _HY_MATRIX_RANDOM_WISHART 03L 60 #define _HY_MATRIX_RANDOM_INVERSE_WISHART 04L 61 #define _HY_MATRIX_RANDOM_MULTINOMIAL 05L 62 63 extern _Trie _HY_MatrixRandomValidPDFs; 64 65 //_____________________________________________________________________________________________ 66 67 class _Formula; 68 /*__________________________________________________________________________________________________________________________________________ */ 69 70 struct _CompiledMatrixData { 71 72 _SimpleFormulaDatum * theStack, 73 * varValues; 74 75 hyFloat * formulaValues; 76 77 long * formulaRefs; 78 bool has_volatile_entries; 79 80 _SimpleList varIndex, 81 formulasToEval; 82 83 }; 84 85 /*__________________________________________________________________________________________________________________________________________ */ 86 87 class _Matrix: public _MathObject { 88 89 // data members 90 public: 91 hyFloat *theData; // matrix elements 92 93 protected: 94 95 // data 96 97 long hDim, vDim, lDim; // matrix physical dimensions; lDim - length of 98 // actual storage allocated 99 100 long* theIndex; // indices of matrix elements in logical storage 101 long* compressedIndex; 102 /** 103 20200821 SLKP to speed sparse caclulations CompressSparseMatrix will create this DENSE index, which has the following structure 104 - Entries in theIndex are expected to be compressed (no -1) and arranged BY ROW 105 - First hDim values: the INDEX of the non-first non zero index for this ROW in theData 106 - Next lDim values: the COLUMN index for the corresponding entry 107 108 [x,1,2,x] 109 [1,2,x,x] 110 [x,x,2,x] 111 [x,x,x,3] 112 113 theIndex : [1,2,4,5,10,15] 114 compressedIndex: [0,2,4,5,1,2,0,1,2,3] 115 116 */ 117 118 private: 119 120 long storageType, // true element matrix(1) or a matrix of pointers(0) which do not need to be deleted 121 // 2, if elements of the matrix are actually formulas which must be initialized to numerics before use 122 // matrices of type two are merely storage tables and can not be operated on directly, i.e their 123 // numerical values are computed first 124 bufferPerRow, // values reflecting internal storage structure for 125 overflowBuffer, // sparse matrices 126 allocationBlock; 127 128 _CompiledMatrixData* 129 cmd; 130 131 HBLObjectRef theValue; // stores evaluated values of the matrix 132 133 static int storageIncrement, // how many percent of full matrix size 134 // to allocate to the matrix storage per increment 135 136 precisionArg, // how many elements in exp series to truncate after 137 138 switchThreshold; // maximum percentage of non-zero elements 139 // to keep the matrix sparse 140 141 static hyFloat truncPrecision; 142 143 // matrix exp truncation precision 144 145 bool _validateCompressedStorage (void) const; 146 147 public: 148 149 150 // constructors 151 152 _Matrix (); // default constructor, doesn't do much 153 154 _Matrix (_String const&, bool, _FormulaParsingContext&, bool use_square_brackets = false); 155 // matrix from a string of the form 156 // {{i11,i12,...}{i21,i22,..}{in1,in2..)}) 157 // or {# rows,<# cols>{i1,j1,expr}{i2,j2,expr}..} 158 // elements may be arbitrary formulas 159 160 _Matrix (long theHDim, long theVDim, bool sparse = false, bool allocateStorage = false); // create an empty matrix of given dimensions; 161 162 163 164 // creates an empty matrix of given dimensions; 165 // the first flag specifies whether it is sparse or not 166 // the second is the storage type -- see below 167 168 // _Matrix (long, long, int, bool); 169 // a test function which generates a random matrix of given dimensions 170 // where the third parameter specifies the percentage of 0 entries and 171 // the first flag indicates how to store the matrix: as spars or usual 172 173 _Matrix ( _Matrix const &); //duplicator 174 175 _Matrix ( _SimpleList const &, long = -1); // make matrix from simple list 176 // the optional argument C (if > 0) tells HyPhy 177 // to make a matrix with C elements per row 178 // if <= 0 - a row matrix is returned 179 180 _Matrix ( _List const &, bool parse_escapes = true); //make string matrix from a list 181 182 _Matrix (const hyFloat *, unsigned long, unsigned long); 183 /* 184 20110616 SLKP 185 added a simple constructor from a list of floating point values 186 first argument: the values (must be at least 2nd arg x 3rd arg long) 187 second argument: how many rows 188 second argument: how many columns 189 190 */ 191 _Matrix (hyFloat, unsigned long, unsigned long); 192 /** 193 20180919 SLKP 194 make an rxc matrix that is constant (each cell is the same) 195 */ 196 197 ~_Matrix (void); //destructor 198 199 virtual void Clear (bool complete = true); //deletes all the entries w/o destroying the matrix 200 virtual void ZeroNumericMatrix (void); //deletes all the entries w/o destroying the matrix 201 202 void Initialize (bool = true); // zeros all matrix structures 203 204 virtual void Serialize (_StringBuffer&,_String&); 205 // write the matrix definition in HBL 206 207 //_____________________________________________________________________________________________ 208 209 210 is_empty(void)211 inline bool is_empty (void) const {return GetVDim () == 0UL || GetHDim () == 0UL;} is_row(void)212 inline bool is_row (void) const { return GetHDim() == 1UL;} is_column(void)213 inline bool is_column (void) const { return GetVDim() == 1UL;} is_square(void)214 inline bool is_square (void) const { return GetVDim() == GetHDim();} is_dense(void)215 inline bool is_dense (void) const {return theIndex == nil;} is_expression_based(void)216 inline bool is_expression_based (void) const {return storageType == _FORMULA_TYPE;} is_numeric(void)217 inline bool is_numeric (void) const {return storageType == _NUMERICAL_TYPE;} is_polynomial(void)218 inline bool is_polynomial (void) const {return storageType == _POLYNOMIAL_TYPE;} has_type(int t)219 inline bool has_type (int t) const {return storageType == t;} 220 221 HBLObjectRef Evaluate (bool replace = true); // evaluates the matrix if contains formulas 222 // if replace is true, overwrites the original 223 224 virtual HBLObjectRef ExecuteSingleOp (long opCode, _List* arguments = nil, _hyExecutionContext* context = _hyDefaultExecutionContext, HBLObjectRef cache = nil); 225 // execute this operation with the list of Args 226 227 HBLObjectRef MAccess (HBLObjectRef, HBLObjectRef, HBLObjectRef cache = nil); 228 // implements the M[i][j] operation for formulas 229 HBLObjectRef MCoord (HBLObjectRef, HBLObjectRef, HBLObjectRef cache = nil); 230 // implements the M[i][j] operation for formulas 231 232 void MStore (long, long, _Formula&, long = -1); 233 void MStore (long, long, HBLObjectRef, long); 234 bool MResolve (HBLObjectRef, HBLObjectRef, long&, long&); 235 // resolve coordiates from two Number arguments 236 237 bool CheckCoordinates ( long&, long&); 238 // validate matrix coordinates 239 240 bool ValidateFormulaEntries (bool (long, long, _Formula*)); 241 // validate matrix coordinates 242 243 void MStore (HBLObjectRef, HBLObjectRef, _Formula&, long = HY_OP_CODE_NONE); 244 // implements the M[i][j]= operation for formulas 245 /* 246 20100811: the last argument provides an op code (-1 = none) 247 to perform on the _Formula argument and the current value in the matrix; 248 this only applies to constant _Formula arguments 249 250 e.g. passing HY_OP_CODE_ADD implements += 251 */ 252 253 void MStore (long, long, HBLObjectRef); 254 void MStore (HBLObjectRef, HBLObjectRef, HBLObjectRef); 255 // implements the M[i][j]= operation for objects 256 virtual HBLObjectRef Compute (void); // returns the numeric value of this matrix 257 258 virtual HBLObjectRef ComputeNumeric (bool = false); // returns the numeric value of this matrix 259 virtual HBLObjectRef RetrieveNumeric (void); // returns the numeric value of this matrix 260 261 virtual void ScanForVariables (_AVLList&, bool inclG = false, _AVLListX* tagger = nil,long weight = 0) const; 262 virtual void ScanForVariables2 (_AVLList&, bool inclG = false, long modelID = -1, bool inclCat = true, _AVLListX* tagger = nil,long weight = 0) const; 263 // scans for all local independent variables on which the matrix depends 264 // and stores them in the list passed as the parameter 265 IsIndependent(void)266 virtual bool IsIndependent (void) { 267 return storageType!=_FORMULA_TYPE; 268 } 269 // used to determine whether the matrix contains references 270 // to other unknowns 271 ObjectClass(void)272 virtual unsigned long ObjectClass (void) const { 273 return MATRIX; 274 } 275 276 _Matrix const& operator = (_Matrix const&); // assignment operation on matrices 277 _Matrix const& operator = (_Matrix const*); // assignment operation on matrices with temp results 278 279 virtual HBLObjectRef Random (HBLObjectRef, HBLObjectRef cache); // reshuffle the matrix 280 281 virtual HBLObjectRef AddObj (HBLObjectRef, HBLObjectRef cache); // addition operation on matrices 282 283 virtual HBLObjectRef SubObj (HBLObjectRef, HBLObjectRef cache); // subtraction operation on matrices 284 285 virtual HBLObjectRef MultObj (HBLObjectRef, HBLObjectRef cache); // multiplication operation on matrices 286 287 virtual HBLObjectRef MultElements (HBLObjectRef, bool elementWiseDivide, HBLObjectRef cache); // element wise multiplication/division operation on matrices 288 289 virtual HBLObjectRef Sum (HBLObjectRef cache); 290 291 _Matrix operator + (_Matrix&); // addition operation on matrices 292 293 _Matrix operator - (_Matrix&); // subtraction operation on matrices 294 295 _Matrix operator * (_Matrix&); // multiplication operation on matrices 296 297 _Matrix operator * (hyFloat); // multiplication of a matrix by a number 298 299 void operator += (_Matrix&); // add/store operation on matrices 300 301 void operator -= (_Matrix&); // subtract/store operation on matrices 302 303 void operator *= (_Matrix&); // multiply/store operation on matrices 304 305 void operator *= (hyFloat); // multiply by a #/store operation on matrices 306 307 void AplusBx (_Matrix&, hyFloat); // A = A + B*x (scalar) 308 309 hyFloat Sqr (hyFloat* _hprestrict_); 310 // square the matrix; takes a scratch vector 311 // of at least lDim doubles 312 // return the maximum absolute element-wise difference between X and X^2 313 314 _List* ComputeRowAndColSums (void); 315 _Matrix* MutualInformation (void); 316 void FillInList (_List&, bool convert_numbers = false) const; 317 // SLKP 20101108: 318 // added a boolean flag to allow numeric matrices 319 // to be implicitly converted to strings 320 321 _Matrix* NeighborJoin (bool, HBLObjectRef cache); 322 _Matrix* MakeTreeFromParent (long, HBLObjectRef cache); 323 _Matrix* ExtractElementsByEnumeration (_SimpleList*,_SimpleList*,bool=false); 324 _Matrix* SimplexSolve (hyFloat = 1.e-6); 325 326 327 // void SqrStrassen (void); 328 // square the matrix by Strassen's Multiplication 329 330 331 _Matrix* Exponentiate (hyFloat scale_to = 0.5, bool check_transition = false, _Matrix * write_here = nil); // exponent of a matrix 332 void Transpose (void); // transpose a matrix 333 _Matrix Gauss (void); // Gaussian Triangularization process 334 HBLObjectRef LUDecompose (void) const; 335 HBLObjectRef CholeskyDecompose (void) const; 336 // added by afyp July 6, 2009 337 HBLObjectRef Eigensystem (HBLObjectRef) const; 338 HBLObjectRef LUSolve (HBLObjectRef) const; 339 HBLObjectRef Inverse (HBLObjectRef cache) const; 340 HBLObjectRef Abs (HBLObjectRef cache); // returns the norm of a matrix 341 // if it is a vector - returns the Euclidean length 342 // otherwise returns the largest element 343 344 hyFloat AbsValue (void) const; 345 346 template <typename CALLBACK> HBLObjectRef ApplyScalarOperation (CALLBACK && functor, HBLObjectRef cache) const; 347 348 // return the matrix of logs of every matrix element 349 350 void SwapRows (const long, const long); 351 long CompareRows (const long, const long); 352 353 hyFloat operator () (long, long) const; // read access to an element in a matrix 354 hyFloat& operator [] (long); // read/write access to an element in a matrix 355 get_dense_numeric_cell(unsigned long r,unsigned long c)356 hyFloat& get_dense_numeric_cell (unsigned long r, unsigned long c) { 357 return theData[r*vDim + c]; 358 } 359 360 void Store (long, long, hyFloat); // write access to an element in a matrix 361 void StoreObject (long, long, _MathObject*, bool dup = false); 362 void StoreObject (long, _MathObject*,bool dup = false); 363 void StoreFormula (long, long, _Formula&, bool = true, bool = true); 364 void NonZeroEntries (_SimpleList&); 365 366 void UpdateDiag (long ,long , _MathObject*); 367 368 void Swap (_Matrix&); // fast swap matrix data 369 friend void SetIncrement (int); // storage parameter access 370 // an auxiliary function which creates an empty 371 // matrix of given dimensions and storage class (normal/sparse) 372 // and storage type (pointer/array) 373 374 friend void DuplicateMatrix (_Matrix*, _Matrix const*); 375 // an auxiliary function which duplicates a matrix 376 377 378 hyFloat MaxElement (char doSum = 0, long * = nil) const; 379 // SLKP 20101022: added an option to return the sum of all elements as an option (doSum = 1) or 380 // the sum of all absolute values (doSum == 2) 381 // returns the largest element's abs value for given matrix 382 // SLKP 20110523: added an option (doSum == 3) to return the largest element (no abs value) 383 // for run modes 0 and 3, if the 2nd argument is non-nil, the index of the winning element will be stored 384 385 hyFloat MinElement (char doAbs = 1, long * = nil); 386 387 // SLKP 20110620: added the 2nd argument to optionally store the index of the smallest element 388 // : added the option to NOT do absolute values 389 // returns the smallest, non-zero element value for given matrix 390 391 bool IsMaxElement (hyFloat); 392 // true if there is an elem abs val of which is greater than the arg 393 // false otherwise 394 395 396 hyFloat MaxRelError(_Matrix&); 397 398 // friend _Matrix IterateStrassen (_Matrix&, _Matrix&); 399 // used in Strassen Squaring 400 401 virtual BaseRef makeDynamic (void) const; // duplicate this object into a dynamic copy 402 403 virtual void Duplicate (BaseRefConst obj); // duplicate this object into a dynamic copy 404 405 virtual BaseRef toStr (unsigned long = 0UL); // convert this matrix to a string 406 407 virtual void toFileStr (FILE*dest, unsigned long = 0UL); 408 409 bool AmISparse (void); 410 411 hyFloat ExpNumberOfSubs (_Matrix*,bool); 412 IsVariable(void)413 virtual bool IsVariable (void) { 414 return storageType != _NUMERICAL_TYPE; 415 } 416 // is this matrix a constant or a variable quantity? 417 418 virtual bool IsConstant (void); 419 IsPrintable(void)420 virtual bool IsPrintable (void) { 421 return storageType != _FORMULA_TYPE; 422 } 423 424 virtual bool Equal (HBLObjectRef); 425 426 void ExportMatrixExp (_Matrix*, FILE*); 427 bool ImportMatrixExp (FILE*); 428 429 hyFloat FisherExact (hyFloat, hyFloat, hyFloat); 430 431 virtual bool HasChanged (bool = false); 432 // have any variables which are referenced by the elements changed? 433 434 virtual unsigned long GetHDim(void)435 GetHDim (void) const{ 436 return hDim; 437 } 438 check_dimension(unsigned long rows,unsigned long columns)439 bool check_dimension (unsigned long rows, unsigned long columns) const { 440 return hDim == rows && vDim == columns; 441 } 442 GetVDim(void)443 unsigned long GetVDim (void) const { 444 return vDim; 445 } GetSize(void)446 unsigned long GetSize (void) const { 447 return lDim; 448 } GetMySize(void)449 long GetMySize (void) { 450 return sizeof(_Matrix)+lDim*(storageType==_NUMERICAL_TYPE?sizeof(hyFloat):sizeof(hyPointer)); 451 } 452 453 void PopulateConstantMatrix (hyFloat); 454 /* SLKP 20090818 455 fill out a numeric matrix with a fixed value 456 if the matrix is sparse, only will out the non-void entries 457 */ 458 459 _Formula* GetFormula (long, long) const; 460 HBLObjectRef GetMatrixCell (long, long, HBLObjectRef cache = nil) const; 461 HBLObjectRef MultByFreqs (long, bool = false); 462 HBLObjectRef EvaluateSimple (_Matrix* existing_receptacle = nil); 463 HBLObjectRef SortMatrixOnColumn (HBLObjectRef, HBLObjectRef cache); 464 HBLObjectRef K_Means (HBLObjectRef, HBLObjectRef cache); 465 HBLObjectRef pFDR (HBLObjectRef, HBLObjectRef cache); // positive false discovery rate 466 HBLObjectRef PoissonLL (HBLObjectRef, HBLObjectRef cache); // log likelihood of a vector of poisson samples given a parameter value 467 468 469 // added by afyp, July 1, 2009 470 HBLObjectRef DirichletDeviate (void); // this matrix used for alpha hyperparameters 471 HBLObjectRef GaussianDeviate (_Matrix &); // " " " " mean hyperparameter, additional argument for variance 472 HBLObjectRef InverseWishartDeviate (_Matrix &); // " " " " rho hyperparameter, additional for phi matrix 473 HBLObjectRef WishartDeviate (_Matrix &, _Matrix &), 474 WishartDeviate (_Matrix &); 475 476 HBLObjectRef MultinomialSample (_Constant*); 477 /* SLKP 20110208: an internal function to draw the multinomial sample 478 479 the matrix _base_ must be 2xN, where each _row_ lists 480 481 value (integer 0 to N-1, but not enforced), probability 482 483 the function will normalize by sum of all the values in the second column 484 485 the constant argument is the number of replicates (M) to draw 486 487 returns an 1xN matrix with counts of how often each value has been drawn 488 489 */ 490 491 492 bool IsValidTransitionMatrix () const; 493 494 bool IsReversible (_Matrix* = nil); 495 // check if the matrix is reversible 496 // if given a base frequencies assumes that rate matrix entries will not be multiplied by freq terms 497 498 bool IsAStringMatrix (void) const; 499 void MakeMeSimple (void); 500 void MakeMeGeneral (void); 501 void ConvertToSimpleList (_SimpleList&); 502 void CompressSparseMatrix (bool, hyFloat*); 503 //prepare the transition probs matrix for exponentiation 504 505 long Hash (long, long) const; // hashing function, which finds matrix 506 // physical element in local storage buffer 507 // returns -1 if insufficient storage 508 // returns a negative number 509 // if element was not found, the number returned 510 // indicates the first available slot 511 fastIndex(void)512 hyFloat* fastIndex(void) const { 513 return (!theIndex)&&(storageType==_NUMERICAL_TYPE)?(hyFloat*)theData:nil; 514 } directIndex(long k)515 inline hyFloat& directIndex(long k) { 516 return theData[k]; 517 } MatrixType(void)518 long MatrixType (void) { 519 return storageType; 520 } 521 ForEach(CALLBACK && cbv,EXTRACTOR && accessor)522 template <typename CALLBACK, typename EXTRACTOR> void ForEach (CALLBACK&& cbv, EXTRACTOR&& accessor) const { 523 if (theIndex) { 524 for (unsigned long i=0UL; i<(unsigned long)lDim; i++) { 525 if (theIndex[i] >= 0L) { 526 cbv (accessor (i), theIndex[i], i); 527 } 528 } 529 } else { 530 for (unsigned long i=0UL; i<(unsigned long)lDim; i++) { 531 cbv (accessor (i), i, i); 532 } 533 } 534 } 535 ForEachCellNumeric(CALLBACK && cbv)536 template <typename CALLBACK> void ForEachCellNumeric (CALLBACK&& cbv) const { 537 if (theIndex) { 538 for (unsigned long i=0UL; i<(unsigned long)lDim; i++) { 539 long idx = theIndex[i]; 540 if (idx >= 0L) { 541 long row = idx / vDim; 542 cbv (theData[i], idx, row, idx - row*vDim); 543 } 544 } 545 } else { 546 for (unsigned long i=0UL, c = 0UL; i<(unsigned long)hDim; i++) { 547 for (unsigned long j=0UL; j<(unsigned long)vDim; j++, c++) { 548 cbv (theData[c], c, i, j); 549 } 550 } 551 } 552 } 553 Any(CALLBACK && cbv,EXTRACTOR && accessor)554 template <typename CALLBACK, typename EXTRACTOR> bool Any (CALLBACK&& cbv, EXTRACTOR&& accessor) const { 555 if (theIndex) { 556 for (unsigned long i=0UL; i<(unsigned long)lDim; i++) { 557 if (theIndex[i] >= 0L) { 558 if (cbv (accessor (i), theIndex[i])) { 559 return true; 560 } 561 } 562 } 563 } else { 564 for (unsigned long i=0UL; i<(unsigned long)lDim; i++) { 565 if (cbv (accessor (i), i)) { 566 return true; 567 } 568 } 569 } 570 return true; 571 } 572 573 bool CheckIfSparseEnough (bool = false, bool copy = true); // check if matrix is sparse enough to justify 574 void Convert2Formulas (void); // converts a numeric matrix to formula-based mx 575 // sparse storage 576 577 void Resize (long); // resize a dense numeric matrix to have more rows 578 579 get_direct(long const index)580 inline hyFloat get_direct (long const index) const { 581 return theData [index]; 582 } 583 get(long const row,long const column)584 inline hyFloat get (long const row, long const column) const { 585 return theData [row * vDim + column]; 586 } 587 set(long const row,long const column)588 inline hyFloat& set (long const row, long const column) { 589 return theData [row * vDim + column]; 590 } 591 592 _String* BranchLengthExpression(_Matrix*, bool); 593 594 void CopyABlock (_Matrix*, long, long, long = 0, long = 0); 595 /* starting at element (row -- 2nd argument, column -- 3rd argument) 596 copy the source matrix (1st argument) row by row 597 598 e.g. if this matrix is 3x4 and the source matrix is 2x2 599 then copying from element 2,2 (0 - based as always) 600 will result in 601 602 [[ x x x x] 603 [ x x x x] 604 [ x x y y]] 605 606 where y is used to denote an element copied from the source 607 4th and 5th arguments override source.hDim and source.vDim, 608 respectively, if they are positive 609 610 Note that both matrices are ASSUMED to be numeric and dense 611 612 NO ERROR CHECKING IS DONE! 613 614 */ 615 /*---------------------------------------------------*/ 616 617 static void CreateMatrix (_Matrix* theMatrix, long theHDim, long theVDim, bool sparse = false, bool allocateStorage = false, bool isFla = false); 618 619 void RecursiveIndexSort (long, long, _SimpleList*); 620 621 622 623 private: 624 625 void internal_to_str (_StringBuffer*, FILE*, unsigned long padding); 626 void SetupSparseMatrixAllocations (void); 627 bool is_square_numeric (bool dense = true) const; 628 629 hyFloat computePFDR (hyFloat, hyFloat); 630 void InitMxVar (_SimpleList& , hyFloat); 631 bool ProcessFormulas (long&, _AVLList&, _SimpleList&, _SimpleList&, _AVLListX&, bool = false, _Matrix* = nil); 632 633 HBLObjectRef PathLogLikelihood (HBLObjectRef, HBLObjectRef cache); 634 /* SLKP: 20100812 635 636 This function assumes that 'this' an 3xK matrix, where each column is of the form 637 A: integer in [0,N-1], B: integer in [0,N-1], T: real >= 0 638 639 and the argument is an NxN RATE matrix for a Markov chain 640 641 The return value is the \sum_ j = 0 ^ {K-1} Prob {A_j -> B_j | T_j} 642 */ 643 644 _Matrix* BranchLengthStencil (void) const; 645 646 //bool IsAStringMatrix (void); 647 void AddMatrix (_Matrix&, _Matrix&, bool sub = false); 648 // aux arithmetic rountines 649 bool AddWithThreshold (_Matrix&, hyFloat); 650 void RowAndColumnMax (hyFloat&, hyFloat&, hyFloat* = nil); 651 void Subtract (_Matrix&, _Matrix&); 652 void Multiply (_Matrix&, hyFloat); 653 void Multiply (_Matrix&, _Matrix const &) const; 654 bool IsNonEmpty (long) const; 655 // checks to see if the i-th position in the storage is non-empty 656 bool CheckDimensions (_Matrix&) const; 657 // compare dims of 2 matrices to see if they can be multiplied 658 long HashBack (long) const; 659 // hashing function, which finds matrix 660 // physical element given local storage 661 void MultbyS (_Matrix&,bool,_Matrix* = nil, hyFloat* = nil); 662 // internal function used in exponentiating sparse matrices 663 664 void Balance (void); // perform matrix balancing; i.e. a norm reduction which preserves the eigenvalues 665 // lifted from balanc function in NR 666 667 void Schur (void); // reduce the matrix to Hessenberg form preserving eigenvalues 668 // lifted from elmhes function in NR 669 670 void EigenDecomp (_Matrix&,_Matrix&) const; // find the eigenvalues of a real matrix 671 // return real and imaginary parts 672 // lifted from hqr in NR 673 674 675 bool AmISparseFast (_Matrix&); 676 677 bool IncreaseStorage (void); 678 // expand the buffer, where elements are held 679 // returns TRUE/FALSE for success/failure 680 681 682 void BreakPoints (long, long, _SimpleList*); 683 void ConvertFormulas2Poly (bool force2numbers=true); 684 void ConvertNumbers2Poly (void); 685 void AgreeObjects (_Matrix&); 686 687 void ClearFormulae (void); // internal reuseable purger 688 void ClearObjects (void); // internal reuseable purger 689 inline 690 _MathObject*& GetMatrixObject(long ind)691 GetMatrixObject (long ind) const { 692 return ((_MathObject**)theData)[ind]; 693 } 694 inline CheckObject(long ind)695 bool CheckObject (long ind) const{ 696 return ((_MathObject**)theData)[ind]!=nil; 697 } 698 699 void SimplexHelper1 (long, _SimpleList&, long, bool, long&, hyFloat&); 700 void SimplexHelper2 (long&, long, hyFloat); 701 void SimplexHelper3 (long,long,long,long); 702 // helper functions for SimplexSolver 703 704 // if nil - matrix stored conventionally 705 706 }; 707 708 /*__________________________________________________________________________________________________________________________________________ */ 709 710 711 /*__________________________________________________________________________________________________________________________________________ */ 712 713 extern _Matrix *GlobalFrequenciesMatrix; 714 // the matrix of frequencies for the trees to be set by block likelihood evaluator 715 extern long ANALYTIC_COMPUTATION_FLAG; 716 717 HBLObjectRef _returnMatrixOrUseCache (long nrow, long ncol, long type, bool is_sparse, HBLObjectRef cache); 718 719 #ifdef _SLKP_USE_AVX_INTRINSICS _avx_sum_4(__m256d const & x)720 inline double _avx_sum_4 (__m256d const & x) { 721 /*__m256d t = _mm256_add_pd (_mm256_shuffle_pd (x, x, 0x0), 722 // (x3,x3,x1,x1) 723 _mm256_shuffle_pd (x, x, 0xf) 724 // (x2,x2,x0,x0); 725 ); 726 return _mm_cvtsd_f64 (_mm_add_pd( 727 _mm256_castpd256_pd128 (t), // (x3+x2,x3+x2) 728 _mm256_extractf128_pd(t,1) // (x1+x0,x0+x1); 729 )); 730 */ 731 __m256d sum = _mm256_hadd_pd(x, x); 732 // sum now has (x[0]+x[1],x[0]+x[1],x[2]+x[3],x[2]+x[3]) 733 return _mm_cvtsd_f64(_mm_add_pd(_mm256_extractf128_pd(sum, 1), _mm256_castpd256_pd128(sum))); 734 /*double __attribute__ ((aligned (32))) array[4]; 735 _mm256_store_pd (array, x); 736 return (array[0]+array[1])+(array[2]+array[3]) ;*/ 737 738 } 739 #endif 740 741 #endif 742