1class Vector {
2   std::vector<double> k;
3   int N;
4public:
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
77class Matrix {
78   int M;
79   std::vector<double*> data;
80public:
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
158Vector 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