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;
}