void ConvertConstraintsToMatrixForm(VectorXi indices, MatrixXd positions, Eigen::SparseMatrix<double> &C, VectorXd &d){
  // Convert the list of fixed indices and their fixed positions to a linear system
  // Hint: The matrix C should contain only one non-zero element per row and d should contain the positions in the correct order.
  C.resize(indices.size()*2, V.rows()*2);
  d.resize(indices.size()*2);
  for(int i=0; i<indices.size(); ++i){
    int idx = indices(i);
    double u = positions(i, 0);
    double v = positions(i, 1); 
    d(i*2+0) = u;
    d(i*2+1) = v;
    C.insert(i*2+0, idx) = 1;
    C.insert(i*2+1, idx) = 1;
  }
}

void computeParameterization(int type){
  assert (1 <= type && type <= 4);
  std::cout << "Computing " << type << "...\n";
  VectorXi fixed_UV_indices;
  MatrixXd fixed_UV_positions;

  SparseMatrix<double> A;
  VectorXd b;
  Eigen::SparseMatrix<double> C;
  VectorXd d;
  // Find the indices of the boundary vertices of the mesh and put them in fixed_UV_indices
  if (!freeBoundary){
    igl::boundary_loop(F, fixed_UV_indices);
    igl::map_vertices_to_circle(V, fixed_UV_indices, fixed_UV_positions);
  }else{
    // FIXME: implement this
  }

  ConvertConstraintsToMatrixForm(fixed_UV_indices, fixed_UV_positions, C, d);

  // Find the linear system for the parameterization (1- Tutte, 2- Harmonic, 3- LSCM, 4- ARAP)
  // and put it in the matrix A.
  // The dimensions of A should be 2#V x 2#V.
  if (type == 1) {
    // Add your code for computing uniform Laplacian for Tutte parameterization
    // Hint: use the adjacency matrix of the mesh
    b = Eigen::VectorXd::Zero(V.rows()*2);
    SparseVector<double> Asum;
    SparseMatrix<double> L, Adj, Adiag;
    igl::adjacency_matrix(F, Adj);
    igl::sum(Adj, 1, Asum);
    igl::diag(Asum, Adiag);
    L = Adj-Adiag;
    // Expand L to [L 0; 0 L] form for A.
    SparseMatrix<double> U, D, Z(L.rows(),L.cols());
    igl::cat(2, L, Z, U);
    igl::cat(2, Z, L, D);
    igl::cat(1, U, D, A);
  }

  if (type == 2) {
    // Add your code for computing cotangent Laplacian for Harmonic parameterization
    // Use can use a function "cotmatrix" from libIGL, but ~~~~***READ THE DOCUMENTATION***~~~~
    b = Eigen::VectorXd::Zero(V.rows()*2);
    igl::cotmatrix(V, F, A);
  }

  if (type == 3) {
    // Add your code for computing the system for LSCM parameterization
    // Note that the libIGL implementation is different than what taught in the tutorial! Do not rely on it!!
  }

  if (type == 4) {
    // Add your code for computing ARAP system and right-hand side
    // Implement a function that computes the local step first
    // Then construct the matrix with the given rotation matrices
  }
  
  assert (A.rows() == V.rows()*2 && A.cols() == V.rows()*2);

  // Solve the linear system.
  // Construct the system as discussed in class and the assignment sheet
  // Use igl::cat to concatenate matrices
  // Use Eigen::SparseLU to solve the system. Refer to tutorial 3 for more detail
  SparseMatrix<double> Lhs, U, D, Ct, Z(C.rows(), C.rows());
  VectorXd rhs(b.size()+d.size()), x;
  Ct = SparseMatrix<double>(C.transpose());
  igl::cat(2, A, Ct, U);
  igl::cat(2, C, Z, D);
  igl::cat(1, U, D, Lhs);
  Lhs.makeCompressed();
  rhs << b, d;
  
  SparseLU<SparseMatrix<double>> solver;
  solver.analyzePattern(Lhs);
  solver.factorize(Lhs);
  solver.compute(Lhs);
  x = solver.solve(rhs);

  // The solver will output a vector
  UV.resize(V.rows(), 2);
  UV.col(0) = x.segment(       0, V.rows());
  UV.col(1) = x.segment(V.rows(), V.rows());
}