1 /* $Id: cuMatrix.cpp 103491 2007-05-04 17:18:18Z kazimird $
2 * ===========================================================================
3 *
4 * PUBLIC DOMAIN NOTICE
5 * National Center for Biotechnology Information
6 *
7 * This software/database is a "United States Government Work" under the
8 * terms of the United States Copyright Act. It was written as part of
9 * the author's official duties as a United States Government employee and
10 * thus cannot be copyrighted. This software/database is freely available
11 * to the public for use. The National Library of Medicine and the U.S.
12 * Government have not placed any restriction on its use or reproduction.
13 *
14 * Although all reasonable efforts have been taken to ensure the accuracy
15 * and reliability of the software and data, the NLM and the U.S.
16 * Government do not and cannot warrant the performance or results that
17 * may be obtained by using this software or data. The NLM and the U.S.
18 * Government disclaim all warranties, express or implied, including
19 * warranties of performance, merchantability or fitness for any particular
20 * purpose.
21 *
22 * Please cite the author in any work or product based on this material.
23 *
24 * ===========================================================================
25 *
26 * Author: David Hurwitz
27 *
28 * File Description:
29 * ported from CDTree app
30 *
31 */
32
33 #include <ncbi_pch.hpp>
34 #include <corelib/ncbi_limits.h>
35 #include <algo/structure/cd_utils/cuMatrix.hpp>
36
37 #include <string.h> // for memset
38
39 BEGIN_NCBI_SCOPE
BEGIN_SCOPE(cd_utils)40 BEGIN_SCOPE(cd_utils)
41
42 void AMatrix_base::MakeArrayBigger(const int RowIndex, const int ColIndex) {
43 //-------------------------------------------------------------
44 // if m_Array is big enough, return true.
45 // if m_Array is too small, make it bigger and return false.
46 //-------------------------------------------------------------
47 // make an array that's big enough (don't let either dimension shrink)
48 int NumRows = Max(m_NumRows, RowIndex+1);
49 int NumCols = Max(m_NumCols, ColIndex+1);
50 AMatrix_base NewMatrix(NumRows, NumCols);
51
52 // copy old array into new array
53 NewMatrix.SlowCopy(*this);
54
55 // free the old array
56 DeAllocate();
57
58 // save new array in place of old array
59 m_Array = NewMatrix.m_Array;
60 m_ColumnFlags = NewMatrix.m_ColumnFlags;
61 m_NumRows = NewMatrix.m_NumRows;
62 m_NumCols = NewMatrix.m_NumCols;
63
64 // pretend new array has been freed
65 // (so that new array isn't freed when leaving this routine)
66 NewMatrix.m_Array = NULL;
67 NewMatrix.m_ColumnFlags = NULL;
68 NewMatrix.m_NumRows = 0;
69 NewMatrix.m_NumCols = 0;
70 }
71
72 // Reduce the size of this, retaining data from the original array
73 // for valid indices of the new array.
Shrink(const int NumRows,const int NumCols)74 bool AMatrix_base::Shrink(const int NumRows, const int NumCols) {
75
76 bool result = false;
77 int i, j;
78
79 if ((NumRows < m_NumRows) || (NumCols < m_NumCols)) {
80 if (NumRows > 0 && NumCols > 0) {
81 //
82 AMatrix_base tmpMatrix(NumRows, NumCols);
83 for (i = 0; i < NumRows; ++i) {
84 for (j = 0; j < NumCols; ++j) {
85 tmpMatrix.m_Array[i][j] = m_Array[i][j];
86 if (i == 0) {
87 tmpMatrix.m_ColumnFlags[j] = m_ColumnFlags[j];
88 }
89 }
90 }
91 DeAllocate();
92 Allocate(NumRows, NumCols);
93 for (i = 0; i < NumRows; ++i) {
94 for (j = 0; j < NumCols; ++j) {
95 m_Array[i][j] = tmpMatrix.m_Array[i][j];
96 if (i == 0) {
97 m_ColumnFlags[j] = tmpMatrix.m_ColumnFlags[j];
98 }
99 }
100 }
101 result = true;
102 }
103 }
104 return result;
105 }
106
Allocate(const int NumRows,const int NumCols)107 void AMatrix_base::Allocate(const int NumRows, const int NumCols) {
108 //--------------------------------------------
109 // allocate memory for array
110 //--------------------------------------------
111 typedef double* DoublePtr;
112
113 // allocate NumRows pointers-to-doubles
114 m_Array = new DoublePtr[NumRows];
115 // for each pointer-to-double
116 for (int i=0; i<NumRows; i++) {
117 // allocate NumCols doubles
118 m_Array[i] = new double[NumCols];
119 }
120
121 // allocate space for column flags, set them to false
122 m_ColumnFlags = new bool[NumCols];
123 memset(m_ColumnFlags, 0, NumCols*sizeof(bool));
124
125 m_NumRows = NumRows;
126 m_NumCols = NumCols;
127 }
128
129
DeAllocate()130 void AMatrix_base::DeAllocate() {
131 //--------------------------------------------
132 // deallocate array memory
133 //--------------------------------------------
134 if (m_Array == NULL) return;
135
136 // for each pointer-to-double
137 for (int i=0; i<m_NumRows; i++) {
138 // free NumCols doubles
139 delete []m_Array[i];
140 m_Array[i] = NULL;
141 }
142 // free NumRows pointer-to-doubles
143 delete []m_Array;
144 m_Array = NULL;
145
146 // free column flags
147 delete []m_ColumnFlags;
148 m_ColumnFlags = NULL;
149
150 m_NumRows = 0;
151 m_NumCols = 0;
152 }
153
154
155
Copy(const AMatrix_base & Matrix)156 void AMatrix_base::Copy(const AMatrix_base& Matrix) {
157 //--------------------------------------------
158 // copy Matrix to this one
159 //--------------------------------------------
160 DeAllocate();
161 Allocate(Matrix.m_NumRows, Matrix.m_NumCols);
162 memcpy(m_Array, Matrix.m_Array, m_NumRows*m_NumCols*sizeof(double));
163 memcpy(m_ColumnFlags, Matrix.m_ColumnFlags, m_NumCols*sizeof(bool));
164 }
165
166
SlowCopy(const AMatrix_base & Matrix)167 void AMatrix_base::SlowCopy(const AMatrix_base& Matrix) {
168 //------------------------------------------------------------------
169 // copy Matrix to this one.
170 // memory for Array must already be allocated.
171 // don't assume this matrix is the same size as Matrix.
172 //------------------------------------------------------------------
173 // if matrices are the same size, go ahead and do a fast copy
174 if ((m_NumRows == Matrix.m_NumRows) && (m_NumCols == Matrix.m_NumCols)) {
175 memcpy(m_Array, Matrix.m_Array, m_NumRows*m_NumCols*sizeof(double));
176 }
177
178 // make sure this matrix is bigger than Matrix
179 assert(m_NumRows >= Matrix.m_NumRows);
180 assert(m_NumCols >= Matrix.m_NumCols);
181
182 // copy Matrix into this matrix
183 for (int i=0; i<Matrix.m_NumRows; i++) {
184 memcpy(m_Array[i], Matrix.m_Array[i], Matrix.m_NumCols*sizeof(double));
185 }
186
187 // copy column flags
188 memcpy(m_ColumnFlags, Matrix.m_ColumnFlags, Matrix.m_NumCols*sizeof(bool));
189 }
190
191
LinearTransform(double b,double m,bool ignoreDiagonal)192 void AMatrix_base::LinearTransform(double b, double m, bool ignoreDiagonal) {
193 int n = GetNumRows();
194 for (int i=0; i<n; ++i) {
195 for (int j=0; j<n; ++j) {
196 if (ignoreDiagonal && i == j) {
197 continue;
198 } else {
199 m_Array[i][j] = m*m_Array[i][j] + b;
200 }
201 }
202 }
203 }
204
205 // Does not assume a symmetric matrix.
GetExtremalEntries(double & max,double & min,bool ignoreDiagonal)206 void AMatrix_base::GetExtremalEntries(double& max, double& min, bool ignoreDiagonal) {
207 int i, j;
208 int n = GetNumRows();
209 double maxEntry=kMin_Double;
210 double minEntry=kMax_Double;
211 for (i=0; i<n; ++i) {
212 for (j=0; j<n; ++j) {
213 if (ignoreDiagonal && i == j) {
214 continue;
215 }
216 if (m_Array[i][j] > maxEntry) {
217 maxEntry = m_Array[i][j];
218 }
219 if (m_Array[i][j] < minEntry) {
220 minEntry = m_Array[i][j];
221 }
222 }
223 }
224 max = maxEntry;
225 min = minEntry;
226 }
227
228 END_SCOPE(cd_utils)
229 END_NCBI_SCOPE
230