| class Vector { |
| std::vector<double> k; |
| int N; |
| public: |
| explicit Vector(int size): k(size) { |
| N = size; |
| for (int i = 0; i < N; i++) |
| k[i] = 0.0; |
| } |
| |
| inline int GetSize() const { return N; } |
| |
| inline double& operator[] (int ix) { |
| CHECK(ix < N && ix >= 0); |
| return k[ix]; |
| } |
| |
| inline const double& operator[] (int ix) const { |
| CHECK(ix < N && ix >= 0); |
| return k[ix]; |
| } |
| |
| inline Vector operator+ (const Vector & other) const { |
| CHECK(other.GetSize() == N); |
| Vector ret(N); |
| for (int i = 0; i < N; i++) |
| ret[i] = k[i] + other[i]; |
| return ret; |
| } |
| |
| inline Vector operator- (const Vector & other) const { |
| CHECK(other.GetSize() == N); |
| Vector ret(N); |
| for (int i = 0; i < N; i++) |
| ret[i] = k[i] - other[i]; |
| return ret; |
| } |
| |
| inline Vector operator* (double factor) const { |
| Vector ret(N); |
| for (int i = 0; i < N; i++) |
| ret[i] = k[i] * factor; |
| return ret; |
| } |
| |
| inline double DotProduct (const Vector & other) const { |
| CHECK(other.GetSize() == N); |
| double ret = 0.0; |
| for (int i = 0; i < N; i++) |
| ret += k[i] * other[i]; |
| return ret; |
| } |
| |
| inline double Len2 () const { |
| double ret = 0.0; |
| for (int i = 0; i < N; i++) |
| ret += k[i] * k[i]; |
| return ret; |
| } |
| |
| inline double Len () const { return sqrt(Len2()); } |
| |
| string ToString () const { |
| char temp[1024] = "("; |
| for (int i = 0; i < N; i++) |
| sprintf(temp, "%s%s%.1lf", temp, i == 0 ? "" : ", ", k[i]); |
| return (std::string)temp + ")"; |
| } |
| |
| inline void Copy(const Vector & from) { |
| CHECK(from.GetSize() == N); |
| for (int i = 0; i < N; i++) |
| k[i] = from[i]; |
| } |
| }; |
| |
| class Matrix { |
| int M; |
| std::vector<double*> data; |
| public: |
| /* |
| * This matrix can grow its "N" i.e. width |
| * See http://en.wikipedia.org/wiki/Matrix_(mathematics) |
| */ |
| |
| Matrix (int M_, int N) { |
| M = M_; |
| for (int i = 0; i < N; i++) { |
| IncN(); |
| } |
| } |
| |
| ~Matrix() { |
| for (int i = 0; i < data.size(); i++) |
| delete [] data[i]; |
| } |
| |
| inline int GetM() const { return M; } |
| |
| inline int GetN() const { return data.size(); } |
| inline void IncN() { |
| double * k = new double[M]; |
| for (int i = 0; i < M; i++) |
| k[i] = 0.0; |
| data.push_back(k); |
| } |
| |
| |
| inline double& At(int i, int j) { |
| CHECK(i < M && i >= 0); |
| CHECK(j < data.size() && j >= 0); |
| // Note the reverse order of indices! |
| return data[j][i]; |
| } |
| |
| inline const double& At(int i, int j) const { |
| CHECK(i < M && i >= 0); |
| CHECK(j < data.size() && j >= 0); |
| // Note the reverse order of indices! |
| return data[j][i]; |
| } |
| |
| Vector MultiplyRight (const Vector & v) const { |
| int N = data.size(); |
| CHECK(v.GetSize() == N); |
| Vector ret(M); |
| for (int i = 0; i < M; i++) |
| for (int j = 0; j < N; j++) |
| ret[i] += v[j] * At(i,j); |
| return ret; |
| } |
| |
| Vector MultiplyLeft (const Vector & v_to_transpose) const { |
| int N = data.size(); |
| CHECK(v_to_transpose.GetSize() == M); |
| Vector ret(N); |
| for (int i = 0; i < M; i++) |
| for (int j = 0; j < N; j++) |
| ret[j] += v_to_transpose[i] * At(i,j); |
| return ret; |
| } |
| |
| string ToString() const { |
| string ret = ""; |
| for (int i = 0; i < M; i++) { |
| ret += "["; |
| for (int j = 0; j < GetN(); j++) { |
| char temp[128] = ""; |
| sprintf(temp, "%s%.1lf", j == 0 ? "" : ", ", At(i,j)); |
| ret += temp; |
| } |
| ret += "]\n"; |
| } |
| return ret; |
| } |
| }; |
| |
| Vector EstimateParameters(const Matrix & perf_m, const Vector & stats_v, double rel_diff, int * iter_count = NULL) |
| { |
| /* |
| * Goal: find Vector<N> parameters: |
| * (perf_m * parameters - stats_v)^2 -> min |
| * With the following limitations: |
| * parameters[i] >= 0 |
| * "Stop rule": |
| * for every i=0...M |
| * -> |stats_v[i] - (perf_m * parameters)[i]| < rel_diff * stats_v[i] |
| * Algorithm used is a gradient descend |
| */ |
| int N = perf_m.GetN(), M = perf_m.GetM(); |
| CHECK(stats_v.GetSize() == M); |
| Vector current(N); |
| |
| // First, let's find out those lines in matrix having only one non-zero coefficient |
| // and if we have any, decrease the dimension of our matrix |
| { |
| int count_easy_param = 0; |
| std::vector<int> parameters_set(N); // N (line#) -> M (row#) of the only nonzero element; -1 if 0/2+ |
| for (int n = 0; n < N; n++) { |
| parameters_set[n] = -1; |
| } |
| // we may detect & assert unresolvable conflicts like the following |
| // |1| |1| |
| // |1| * p = |2| |
| // |1| |3| |
| |
| // Find out those 0000*0000 lines |
| for (int m = 0; m < M; m++) { |
| int zero_id = -1; // -1 = not found, -2 = found twice, else - idx |
| for (int n = 0; n < N; n++) { |
| const double & m_n = perf_m.At(m,n); |
| if (m_n != 0.0) { |
| if (zero_id == -1) { |
| zero_id = n; |
| } else if (zero_id >= 0) { |
| zero_id = -2; |
| break; |
| } |
| } |
| } |
| if (zero_id >= 0) { |
| CHECK(parameters_set[zero_id] == -1); |
| parameters_set[zero_id] = m; |
| count_easy_param++; |
| current[zero_id] = stats_v[m] / perf_m.At(m, zero_id); |
| } |
| } |
| |
| if (count_easy_param > 0) { |
| //printf("tmp = %s\n", current.ToString().c_str()); |
| |
| //printf("\n\nDecreasing the dimensions!\n"); |
| std::vector<int> new_m_to_old(M - count_easy_param), |
| new_n_to_old(N - count_easy_param); |
| for (int m = 0; m < M - count_easy_param; m++) { |
| // see increments later |
| new_m_to_old[m] = m; |
| } |
| int cur_n = 0; |
| for (int n = 0; n < N; n++) { |
| if (parameters_set[n] == -1) { |
| new_n_to_old[cur_n] = n; |
| cur_n++; |
| } else { |
| for (int m = parameters_set[n]; m < M - count_easy_param; m++) { |
| new_m_to_old[m]++; |
| } |
| } |
| } |
| |
| Vector auto_stats = perf_m.MultiplyRight(current), // we have these stats for sure ... |
| diff_stats = stats_v - auto_stats; // ... and we need to get these stats more |
| |
| // Formulate the sub-problem (sub-matrix & sub-stats) |
| Matrix new_m(M - count_easy_param, N - count_easy_param); |
| Vector new_stats(M - count_easy_param); |
| for (int m = 0; m < M - count_easy_param; m++) { |
| new_stats[m] = diff_stats[new_m_to_old[m]]; |
| CHECK(new_stats[m] >= 0.0); |
| for (int n = 0; n < N - count_easy_param; n++) |
| new_m.At(m,n) = perf_m.At(new_m_to_old[m], new_n_to_old[n]); |
| } |
| |
| // Solve the submatrix and put the results back |
| Vector new_param = EstimateParameters(new_m, new_stats, rel_diff, iter_count); |
| for (int n = 0; n < N - count_easy_param; n++) { |
| current[new_n_to_old[n]] = new_param[n]; |
| } |
| return current; |
| } |
| } |
| |
| bool stop_condition = false; |
| double prev_distance = stats_v.Len2(); |
| |
| //printf("Initial distance = %.1lf\n", prev_distance); |
| while (stop_condition == false) { |
| (*iter_count)++; |
| Vector m_by_curr(M); |
| m_by_curr.Copy(perf_m.MultiplyRight(current)); |
| Vector derivative(N); // this is actually "-derivative/2" :-) |
| derivative.Copy(perf_m.MultiplyLeft(stats_v - m_by_curr)); |
| |
| if (derivative.Len() > 1000.0) { |
| derivative = derivative * (1000.0 / derivative.Len()); |
| } |
| |
| // Descend according to the gradient |
| { |
| Vector new_curr(N); |
| double step = 1024.0; |
| double new_distance; |
| for (int i = 0; i < N; i++) { |
| if (current[i] + derivative[i] * step < 0.0) { |
| step = - current[i] / derivative[i]; |
| } |
| } |
| do { |
| new_curr = current + derivative*step; |
| step /= 2.0; |
| new_distance = (perf_m.MultiplyRight(new_curr) - stats_v).Len2(); |
| if (step < 0.00001) { |
| stop_condition = true; |
| } |
| (*iter_count)++; |
| } while (new_distance >= prev_distance && !stop_condition); |
| |
| prev_distance = new_distance; |
| current = new_curr; |
| //printf ("Dist=%.1lf, cur=%s\n", new_distance, current.ToString().c_str()); |
| } |
| |
| // Check stop condition |
| if (stop_condition == false) |
| { |
| stop_condition = true; |
| Vector stats_est(M); |
| stats_est.Copy(perf_m.MultiplyRight(current)); |
| for (int i = 0; i < M; i++) |
| if (fabs(stats_v[i] - stats_est[i]) > rel_diff * stats_v[i]) |
| stop_condition = false; |
| } |
| } |
| |
| return current; |
| } |