Project a vector of reals onto a simplex in julia

projectsimplex.jl

# Figure 1 algorithm from http://icml2008.cs.helsinki.fi/papers/361.pdf
# project a real vector onto a simplex in nlogn time
# v is the vector to project
# z is the value to which the projected vector must sum, one by default
function projectsimplex{T <: Real}(v::Array{T, 1}, z::T=one(T))
	n=length(v)
	µ = sort(v, rev=true)
	#finding ρ could be improved to avoid so much temp memory allocation
	ρ = maximum((1:n)[µ - (cumsum(µ) .- z) ./ (1:n) .> zero(T)]) 
	θ = (sum(µ[1:ρ]) - z)/ρ
	max(v .- zero(T), zero(T))
end

i suppose the last line should read

max(v .- θ, zero(T))