1 class Vector { 2 std::vector<double> k; 3 int N; 4 public: 5 explicit Vector(int size): k(size) { 6 N = size; 7 for (int i = 0; i < N; i++) 8 k[i] = 0.0; 9 } 10 11 inline int GetSize() const { return N; } 12 13 inline double& operator[] (int ix) { 14 CHECK(ix < N && ix >= 0); 15 return k[ix]; 16 } 17 18 inline const double& operator[] (int ix) const { 19 CHECK(ix < N && ix >= 0); 20 return k[ix]; 21 } 22 23 inline Vector operator+ (const Vector & other) const { 24 CHECK(other.GetSize() == N); 25 Vector ret(N); 26 for (int i = 0; i < N; i++) 27 ret[i] = k[i] + other[i]; 28 return ret; 29 } 30 31 inline Vector operator- (const Vector & other) const { 32 CHECK(other.GetSize() == N); 33 Vector ret(N); 34 for (int i = 0; i < N; i++) 35 ret[i] = k[i] - other[i]; 36 return ret; 37 } 38 39 inline Vector operator* (double factor) const { 40 Vector ret(N); 41 for (int i = 0; i < N; i++) 42 ret[i] = k[i] * factor; 43 return ret; 44 } 45 46 inline double DotProduct (const Vector & other) const { 47 CHECK(other.GetSize() == N); 48 double ret = 0.0; 49 for (int i = 0; i < N; i++) 50 ret += k[i] * other[i]; 51 return ret; 52 } 53 54 inline double Len2 () const { 55 double ret = 0.0; 56 for (int i = 0; i < N; i++) 57 ret += k[i] * k[i]; 58 return ret; 59 } 60 61 inline double Len () const { return sqrt(Len2()); } 62 63 string ToString () const { 64 char temp[1024] = "("; 65 for (int i = 0; i < N; i++) 66 sprintf(temp, "%s%s%.1lf", temp, i == 0 ? "" : ", ", k[i]); 67 return (std::string)temp + ")"; 68 } 69 70 inline void Copy(const Vector & from) { 71 CHECK(from.GetSize() == N); 72 for (int i = 0; i < N; i++) 73 k[i] = from[i]; 74 } 75 }; 76 77 class Matrix { 78 int M; 79 std::vector<double*> data; 80 public: 81 /* 82 * This matrix can grow its "N" i.e. width 83 * See http://en.wikipedia.org/wiki/Matrix_(mathematics) 84 */ 85 86 Matrix (int M_, int N) { 87 M = M_; 88 for (int i = 0; i < N; i++) { 89 IncN(); 90 } 91 } 92 93 ~Matrix() { 94 for (int i = 0; i < data.size(); i++) 95 delete [] data[i]; 96 } 97 98 inline int GetM() const { return M; } 99 100 inline int GetN() const { return data.size(); } 101 inline void IncN() { 102 double * k = new double[M]; 103 for (int i = 0; i < M; i++) 104 k[i] = 0.0; 105 data.push_back(k); 106 } 107 108 109 inline double& At(int i, int j) { 110 CHECK(i < M && i >= 0); 111 CHECK(j < data.size() && j >= 0); 112 // Note the reverse order of indices! 113 return data[j][i]; 114 } 115 116 inline const double& At(int i, int j) const { 117 CHECK(i < M && i >= 0); 118 CHECK(j < data.size() && j >= 0); 119 // Note the reverse order of indices! 120 return data[j][i]; 121 } 122 123 Vector MultiplyRight (const Vector & v) const { 124 int N = data.size(); 125 CHECK(v.GetSize() == N); 126 Vector ret(M); 127 for (int i = 0; i < M; i++) 128 for (int j = 0; j < N; j++) 129 ret[i] += v[j] * At(i,j); 130 return ret; 131 } 132 133 Vector MultiplyLeft (const Vector & v_to_transpose) const { 134 int N = data.size(); 135 CHECK(v_to_transpose.GetSize() == M); 136 Vector ret(N); 137 for (int i = 0; i < M; i++) 138 for (int j = 0; j < N; j++) 139 ret[j] += v_to_transpose[i] * At(i,j); 140 return ret; 141 } 142 143 string ToString() const { 144 string ret = ""; 145 for (int i = 0; i < M; i++) { 146 ret += "["; 147 for (int j = 0; j < GetN(); j++) { 148 char temp[128] = ""; 149 sprintf(temp, "%s%.1lf", j == 0 ? "" : ", ", At(i,j)); 150 ret += temp; 151 } 152 ret += "]\n"; 153 } 154 return ret; 155 } 156 }; 157 158 Vector EstimateParameters(const Matrix & perf_m, const Vector & stats_v, double rel_diff, int * iter_count = NULL) 159 { 160 /* 161 * Goal: find Vector<N> parameters: 162 * (perf_m * parameters - stats_v)^2 -> min 163 * With the following limitations: 164 * parameters[i] >= 0 165 * "Stop rule": 166 * for every i=0...M 167 * -> |stats_v[i] - (perf_m * parameters)[i]| < rel_diff * stats_v[i] 168 * Algorithm used is a gradient descend 169 */ 170 int N = perf_m.GetN(), M = perf_m.GetM(); 171 CHECK(stats_v.GetSize() == M); 172 Vector current(N); 173 174 // First, let's find out those lines in matrix having only one non-zero coefficient 175 // and if we have any, decrease the dimension of our matrix 176 { 177 int count_easy_param = 0; 178 std::vector<int> parameters_set(N); // N (line#) -> M (row#) of the only nonzero element; -1 if 0/2+ 179 for (int n = 0; n < N; n++) { 180 parameters_set[n] = -1; 181 } 182 // we may detect & assert unresolvable conflicts like the following 183 // |1| |1| 184 // |1| * p = |2| 185 // |1| |3| 186 187 // Find out those 0000*0000 lines 188 for (int m = 0; m < M; m++) { 189 int zero_id = -1; // -1 = not found, -2 = found twice, else - idx 190 for (int n = 0; n < N; n++) { 191 const double & m_n = perf_m.At(m,n); 192 if (m_n != 0.0) { 193 if (zero_id == -1) { 194 zero_id = n; 195 } else if (zero_id >= 0) { 196 zero_id = -2; 197 break; 198 } 199 } 200 } 201 if (zero_id >= 0) { 202 CHECK(parameters_set[zero_id] == -1); 203 parameters_set[zero_id] = m; 204 count_easy_param++; 205 current[zero_id] = stats_v[m] / perf_m.At(m, zero_id); 206 } 207 } 208 209 if (count_easy_param > 0) { 210 //printf("tmp = %s\n", current.ToString().c_str()); 211 212 //printf("\n\nDecreasing the dimensions!\n"); 213 std::vector<int> new_m_to_old(M - count_easy_param), 214 new_n_to_old(N - count_easy_param); 215 for (int m = 0; m < M - count_easy_param; m++) { 216 // see increments later 217 new_m_to_old[m] = m; 218 } 219 int cur_n = 0; 220 for (int n = 0; n < N; n++) { 221 if (parameters_set[n] == -1) { 222 new_n_to_old[cur_n] = n; 223 cur_n++; 224 } else { 225 for (int m = parameters_set[n]; m < M - count_easy_param; m++) { 226 new_m_to_old[m]++; 227 } 228 } 229 } 230 231 Vector auto_stats = perf_m.MultiplyRight(current), // we have these stats for sure ... 232 diff_stats = stats_v - auto_stats; // ... and we need to get these stats more 233 234 // Formulate the sub-problem (sub-matrix & sub-stats) 235 Matrix new_m(M - count_easy_param, N - count_easy_param); 236 Vector new_stats(M - count_easy_param); 237 for (int m = 0; m < M - count_easy_param; m++) { 238 new_stats[m] = diff_stats[new_m_to_old[m]]; 239 CHECK(new_stats[m] >= 0.0); 240 for (int n = 0; n < N - count_easy_param; n++) 241 new_m.At(m,n) = perf_m.At(new_m_to_old[m], new_n_to_old[n]); 242 } 243 244 // Solve the submatrix and put the results back 245 Vector new_param = EstimateParameters(new_m, new_stats, rel_diff, iter_count); 246 for (int n = 0; n < N - count_easy_param; n++) { 247 current[new_n_to_old[n]] = new_param[n]; 248 } 249 return current; 250 } 251 } 252 253 bool stop_condition = false; 254 double prev_distance = stats_v.Len2(); 255 256 //printf("Initial distance = %.1lf\n", prev_distance); 257 while (stop_condition == false) { 258 (*iter_count)++; 259 Vector m_by_curr(M); 260 m_by_curr.Copy(perf_m.MultiplyRight(current)); 261 Vector derivative(N); // this is actually "-derivative/2" :-) 262 derivative.Copy(perf_m.MultiplyLeft(stats_v - m_by_curr)); 263 264 if (derivative.Len() > 1000.0) { 265 derivative = derivative * (1000.0 / derivative.Len()); 266 } 267 268 // Descend according to the gradient 269 { 270 Vector new_curr(N); 271 double step = 1024.0; 272 double new_distance; 273 for (int i = 0; i < N; i++) { 274 if (current[i] + derivative[i] * step < 0.0) { 275 step = - current[i] / derivative[i]; 276 } 277 } 278 do { 279 new_curr = current + derivative*step; 280 step /= 2.0; 281 new_distance = (perf_m.MultiplyRight(new_curr) - stats_v).Len2(); 282 if (step < 0.00001) { 283 stop_condition = true; 284 } 285 (*iter_count)++; 286 } while (new_distance >= prev_distance && !stop_condition); 287 288 prev_distance = new_distance; 289 current = new_curr; 290 //printf ("Dist=%.1lf, cur=%s\n", new_distance, current.ToString().c_str()); 291 } 292 293 // Check stop condition 294 if (stop_condition == false) 295 { 296 stop_condition = true; 297 Vector stats_est(M); 298 stats_est.Copy(perf_m.MultiplyRight(current)); 299 for (int i = 0; i < M; i++) 300 if (fabs(stats_v[i] - stats_est[i]) > rel_diff * stats_v[i]) 301 stop_condition = false; 302 } 303 } 304 305 return current; 306 } 307
