function w = budget_voting(X)
  # usage: budget_voting(X)
  #
  # X is a matrix with voters as rows and candidates as columns; each entry
  # reflects a relative budget assigned by the voter to the candidate.
  #
  # The function returns a column vector of relative budget assignments for
  # each candidate such that the sum of Manhattan distances is minimized
  # and the total budget is spent. The Euclidean distance and, as a
  # fallback, the 1-dimensional median are used as tie-breaker.

  # Xt is X transposed, i.e. candidates as rows and voters as columns.
  # Xtn is Xn normalized such that each voter assigns a total budget of 1.
  # n is the number of voters and d is the number of candidates.
  # avg is the arithmetic mean of all ballots.
  Xt = transpose(X);
  Xtn = Xt ./ sum(Xt, 1);
  n = columns(Xtn); d = rows(Xtn);
  avg = mean(Xtn, 2);

  # Handle the 1-dimensional case:
  if rank(Xtn - avg) <= 1; w = median(Xtn, 2); return; endif

  # lb are the lower and ub are the upper bounds.
  lb = zeros(d, 1); ub = ones(d, 1);

  # cnd1(x) is zero if sum(x)=1.
  # cnd1_grad_t is gradient of cnd1 as a row vector.
  cnd1 = @(x) sum(x, 1) - 1;
  cnd1_grad_t = @(x) ones(1, d);

  # obj1 is the first objective function to minimize.
  # obj1_grad is the gradient of obj1 as a column vector.
  obj1 = @(x) sum(sum(abs(x-Xtn), 1), 2);
  obj1_grad = @(x) sum(sign(x-Xtn), 2);

  # Optimize for obj1:
  w = sqp(
    avg, {obj1, obj1_grad},  # start value, objective function, gradient
    {cnd1, cnd1_grad_t},     # equality constraints
    [], lb, ub,              # inequality constraints, lower/upper bounds
    1000                     # max iterations
  );

  # Normalize in case of numerical errors and determine optimal error:
  w = w / sum(w, 1);
  opt1 = obj1(w);

  # cnd2(x) >= 0 for all rows if obj1(x) is minimal (within tolerance tol).
  # cnd2_grad_m is gradient of cnd2 as a matrix.
  tol = n * d * eps;
  cnd2 = @(x) tol + opt1 - obj1(x);
  cnd2_grad_m = @(x) -obj1_grad(x)';

  # obj2 is the second objective function to minimize.
  # obj2_grad is the gradient of obj2 as a column vector.
  obj2 = @(x) sum(sqrt(sum((x-Xtn).^2, 1)), 2);
  obj2_grad = @(x) sum((x-Xtn) ./ sqrt(sum((x-Xtn).^2, 1)), 2);

  # Optimize for obj2:
  w = sqp(
    w, {obj2, obj2_grad},  # start value, objective function, gradient
    {cnd1, cnd1_grad_t},   # equality constraints
    {cnd2, cnd2_grad_m},   # inequality constraints
    lb, ub,                # lower/upper bounds
    1000                   # max iterations
  );

  # Normalize in case of numerical errors:
  w = w / sum(w, 1);
end
