Don't wrap round

[dpsanders]

This seems like a gotcha (indeed it caught us out when manipulating images):

julia> x = N0f8(0.5)
0.502N0f8

julia> x + x
0.0N0f8

[timholy]

#41. FWIW this happens in most image-processing frameworks, including Matlab (which uses saturating arithmetic) and Python (which uses overflow arithmetic like us). That said, I am tempted to blaze a new trail to make this better.

There was an epic discussion in #143. We seem close to re-addressing this shortly, with the following proposals:

• No supported arithmetic, convert immediately to float
• Overflow (wrap-around) arithmetic like it is now
• Saturating arithmetic like Matlab
• Checked arithmetic (you'll get an error from that operation)
• Widening arithmetic like in WIP/RFC/Scary: allow signed Normed numbers #143 (aside from supporting signed arithmetic, the key element of that proposal was a 2-fold widening, but no further widening from that point forward---the types are marked as "fully widened" and henceforth use overflow arithmetic)

Each proposal has its strengths and weaknesses. I'd be curious to hear your thoughts. I know @kimikage is planning to solicit input on discourse, but since you posted it...

[dpsanders]

Thanks. I guess we were expecting it to saturate, but I see that it's a much harder (or insoluble) problem than it initially seems.

One option would be to implement +(:saturate, x , y) or +(Saturate(), x, y) etc., but of course this is pretty cumbersome. (And I guess have plain + throw and error and inform you about these options.)

[timholy]

It's not impossible, it's just not obvious it's a good idea. Stefan has posted on the evils of saturating artihmetic. (E.g., x + x - x != x for many x.)

[timholy]

Also much of it is already implemented (by @kimikage), see

FixedPointNumbers.jl/src/FixedPointNumbers.jl

Lines 218 to 287 in 595f7a7

	# saturating arithmetic
	saturating_neg(x::X) where {X <: FixedPoint} = X(~min(x.i - true, x.i), 0)
	saturating_neg(x::X) where {X <: FixedPoint{<:Unsigned}} = zero(X)
	saturating_abs(x::X) where {X <: FixedPoint} =
	X(ifelse(signbit(abs(x.i)), typemax(x.i), abs(x.i)), 0)
	saturating_add(x::X, y::X) where {X <: FixedPoint} =
	X(x.i + ifelse(x.i < 0, max(y.i, typemin(x.i) - x.i), min(y.i, typemax(x.i) - x.i)), 0)
	saturating_add(x::X, y::X) where {X <: FixedPoint{<:Unsigned}} = X(x.i + min(~x.i, y.i), 0)
	saturating_sub(x::X, y::X) where {X <: FixedPoint} =
	X(x.i - ifelse(x.i < 0, min(y.i, x.i - typemin(x.i)), max(y.i, x.i - typemax(x.i))), 0)
	saturating_sub(x::X, y::X) where {X <: FixedPoint{<:Unsigned}} = X(x.i - min(x.i, y.i), 0)
	saturating_mul(x::X, y::X) where {X <: FixedPoint} = clamp(float(x) * float(y), X)
	saturating_fdiv(x::X, y::X) where {X <: FixedPoint} =
	clamp(floattype(X)(x.i) / floattype(X)(y.i), X)
	# checked arithmetic
	checked_neg(x::X) where {X <: FixedPoint} = checked_sub(zero(X), x)
	function checked_abs(x::X) where {X <: FixedPoint}
	abs(x.i) >= 0 || throw_overflowerror_abs(x)
	X(abs(x.i), 0)
	end
	function checked_add(x::X, y::X) where {X <: FixedPoint}
	r, f = Base.Checked.add_with_overflow(x.i, y.i)
	z = X(r, 0) # store first
	f && throw_overflowerror(:+, x, y)
	z
	end
	function checked_sub(x::X, y::X) where {X <: FixedPoint}
	r, f = Base.Checked.sub_with_overflow(x.i, y.i)
	z = X(r, 0) # store first
	f && throw_overflowerror(:-, x, y)
	z
	end
	function checked_mul(x::X, y::X) where {X <: FixedPoint}
	z = float(x) * float(y)
	typemin(X) - eps(X)/2 <= z < typemax(X) + eps(X)/2 || throw_overflowerror(:*, x, y)
	z % X
	end
	function checked_fdiv(x::X, y::X) where {T, X <: FixedPoint{T}}
	y === zero(X) && throw(DivideError())
	z = floattype(X)(x.i) / floattype(X)(y.i)
	if T <: Unsigned
	z < typemax(X) + eps(X)/2 || throw_overflowerror(:/, x, y)
	else
	typemin(X) - eps(X)/2 <= z < typemax(X) + eps(X)/2 || throw_overflowerror(:/, x, y)
	end
	z % X
	end
	# default arithmetic
	const DEFAULT_ARITHMETIC = :wrapping
	for (op, name) in ((:-, :neg), (:abs, :abs))
	f = Symbol(DEFAULT_ARITHMETIC, :_, name)
	@eval begin
	$op(x::X) where {X <: FixedPoint} = $f(x)
	end
	end
	for (op, name) in ((:+, :add), (:-, :sub), (:*, :mul))
	f = Symbol(DEFAULT_ARITHMETIC, :_, name)
	@eval begin
	$op(x::X, y::X) where {X <: FixedPoint} = $f(x, y)
	end
	end
	/(x::X, y::X) where {X <: FixedPoint} = checked_fdiv(x, y) # force checked arithmetic

[kimikage]

We should eventually decide on a final course of action on the overflow issue, but we haven't made any significant progress in the last six months at least. 😭

Apart from the overflow issues, I was unhappy with the implementation of arithmetic (e.g. ::Normed * ::Normed). So I decided to support checked, wrapping and saturating arithmetic as part of the improvement of arithmetic implementations.

After I complete the implementation of checked arithmetic (cf. #41 (comment)), I'm going to open a thread in discourse to decide which is the default arithmetic (on a per-method basis) in v0.9.0.

That decision in v0.9.0 will only be the first step and should not be the final goal. However, I believe it is useful to have real working codes in order to determine the goal. Benchmarking is extremely important because overflow solutions are closely related to speed and memory performance.

[kimikage]

Incidentally, I "personally" disagree with making saturating arithmetic the default. However, I think there is a rationale for using saturation arithmetic for Gray{N0f8} and RGB{N0f8}, not N0f8.

I'm also preparing to experiment with the three types of arithmetic for ColorVectorSpace. https://github.com/kimikage/CheckedColorArithmetic.jl

[kimikage]

@timholy, I would appreciate it if you could reply to JuliaMath/CheckedArithmetic.jl#6 when you have a time.

[timholy]

As you think about the options, let me make a pitch again for the type of solution in #143, since you haven't mentioned it as something you're actively considering and because I do think it's one of the more attractive options.

The naive motivation is the following: if UInt8 + UInt8 might overflow, why not just return a UInt64? Well, that's a lot bigger and eats up a lot more memory, a serious negative if you're processing large images. OK, so how about a UInt16? The problem there is that when you then go to add a UInt16 to anything else, do you return a UInt32? That quickly gets out of hand---it makes arithmetic essentially non-inferrable.

So the core idea in #143 is to split into two worlds: "representation types" and "computed types." "Representation types" are what might be created, say, by reading an image from a file---they are used for your typical "starting image." Computed types are produced as the result of any math operation (well, any operation that tries to stay in the world of Fixed* numbers). Continuing the UInt analogy, instead of UInt8 + UInt8 -> UInt16, we instead have UInt8 + UInt8 -> WUInt16, where W marks it as a computed type. The key proposal is that representation types widen to computed types, but that computed types are "terminal" and used wrap-around (overflow) arithmetic. What that means is that if you're playing around at the command line and add two images and then send the result to a viewer, you should get exactly what you expect; if we make the computed type signed, then the same is true of the difference between two images. That latter is especially important: in image processing, we look at the difference image all the time, and making that both easy and accurate seems like an important goal.

The biggest strength of the proposal is that naive operations "just work," which is true of only one of the other proposals, having all math operations immediately convert to float. Unlike the latter, it also keeps resulting values relatively compact. The weaknesses of the proposal are that it effectively doubles the number of types we have to support, and it could lull people into false complacency: my proposed strategy is great for interactive exploration where you compute the sum or difference of a pair of images, but if you sum 1000 UInt8-represented images you need to allocate the result using a safer type.

Why don't other languages/frameworks do it this way? Because they don't have the types---they are limited to machine numbers, and so they can't distinguish a terminally-widened WUInt16 from a representational UInt16. But we don't have such limitations.

On balance I think having things "just work" for interactive exploration, while not compromising performance and accuracy, is probably the most important goal. Hence my two favorite options (in the abstract, without any practical experience) are "all math immediately converts to float" and "limited doubling of size" (aka, #143). I'm less fond of the saturating, "simple" wrap-around, or checked arithmetic options, which are the main ones you have been promoting. Hence I wanted to remind you of these options.

[johnnychen94]

Unsure how much effort is needed to achieve this, from the user's side I like the idea of splitting it apart to representation types and computation types. If I get the idea correctly, the representation types are mostly exclusive to IO package developers and are transparent to normal users.

Haven't taken a time to experiment with the new checked-math functions developed by @kimikage, though.

[timholy]

Yes, the idea is that we want to be able to represent any "input" image without conversion; reinterpret is fine, I'm talking about the bit-level description here. This is particularly important if your image is a 1TB file that you import via memory-mapping: you must be able to represent the meaning of that file using its native bit layout. But once you're working with it in memory you're free to choose different types for the result of math operations.

[kimikage]

The fundamental problem is that what is a "representation type" and what is a "computed type" depend on the use case. For me, N0f8 is sometimes a representation type and sometimes a computed type.

If you think it's an exception, please consider the case of ::N0f8 * ::N0f8. This can cause a loss of accuracy, but not an overflow. If you convert N0f8 to a wider type, you lose the optimization opportunity described above.

In fact, there must be both cases where you want to widen and cases where you don't want to widen. The two requirements are completely antithetical and there is no universal solution. The "neutral" approach is to let the user specify the way. We can add widening_* methods, but users can still convert the types manually.

[timholy]

Sure, but in cases where widening isn't necessary then you don't have to do it. Maybe instead of calling them "representation" and "computed" we could call them "plain" and "marked." "plain" types can widen, if the operation requires it to avoid overflow; "marked" types cannot.

[kimikage]

You seem to be addressing only a single operation. I'm illustrating a case like (a::N0f8 - b::N0f8) * c::N0f8. How do I clear the "marked"?

Edit:
Of course, this is about regressions, so I'm assuming that a - b is working correctly (even if it's wrapped-around). (1 - b) * c is actually used in interpolation etc.

Edit2:

Edit: you can always use % N0f8 for interpolation, or better just write that as complement(b).

Yes, I know the complement. However, it doesn't seem to me that calling additional functions intentionally is much different than leaving it up to the user to write float(x) or N24f8(x). In some cases, just widening 8-bit types into 16-bit types doesn't "just work".

[timholy]

You don't. The answer might be negative. But in any case, once marked, always marked.

Edit: you can always use % N0f8 for interpolation, or better just write that as complement(b).

[kimikage]

The other major problem is that the 16-bit types are too rich for the 8-bit types. There is little benefit in using FixedPoint if it is allowed to be extended to 16-bit types. I've demonstrated before that Float16 is also never too slow if we give up full compliance with IEEE 754. If you are still concerned about speed (that is already a contradiction when using the 16-bit type, though), then the JuliaImages community should think more seriously about supporting BFloat16.

Also, in image processing, there are 8 bits that are not used in RGB24. Even if you add 2 bits to each channel, there is still 2 bits left. Both NaN and Inf can be represented if they are not channel-wise. Of course, the type is overly application-specific, but I think the types you're trying to add are overly application-specific, too.

[timholy]

I'm fine with a fast Float16; I basically consider #143 and #147 different approaches to the same thing. BFloat16s is pretty interesting too. I guess there we can do some manipulations in Float32 and then strip bits?

Not all images are RGB24. My "real" work is mostly grayscale.

[timholy]

Note: we don't necessarily have to choose, after all we have 4 (way more, actually) untapped binary addition operators ⨣ ⨤ ⨥ ⨦ 😝

[kimikage]

BTW, I suggested in #167 to use f as a marker, rather than introducing new marked types. In any case, though, we will need to add signed Normed. (I don't think signed Normed is application-specific.)

[nlw0]

I was considering to open a new issue, then I saw this one that seems related and fresh.

I'm working with N0f8 images, and I have a windowing function which happens to be N1f7. Ideally when I do the operation I was expecting the result to be N0f8 * N1f7 = N1f15. What we get instead is N8f8. I imagine the rationale is just to keep the same largest accuracy, what makes N0f8 * N0f8 return N0f8?

I'm myself keen on preserving accuracy, with the best speed. Expanding from 8 bits to 16 bits is as far as I know the natural result in the low level, and in my example above it becomes a simple multiplication, while ending in N8f8 actually should involve a shift-right, if I'm not mistaken.

What is the rationale for these conversions right now, and is there any plans to support a type that preserves accuracy?

[timholy]

The simplest rationale is that this is what Base does:

julia> x = 0x02
0x02

julia> x*x
0x04

Regarding your proposal to return N1f15, see the "getting out of hand" remark in #227 (comment). I'd say you should do float(x) * float(y) if you're keen to preserve accuracy.

[kimikage]

Another important point is that the conversion of Normed's f is generally lossy. Of course, the current promotion rules are not perfect, but they are not so bad either.

BTW, the Normed-->Normed conversions have room for improvement. (cf. #140)
(Edit: As #140 shows, N2f14 is much better than N1f15. 😄)

Edit:
I am also prepared to support widemul if needed. However, it is something like widemul(::N0f8, ::N0f8)::N8f8. I don't want to implement widemul for types with different fs. 😅

[nlw0]

I understand for most applications preserving a constant precision is the way to go, especially if you have long expressions and it would "get out of hand". This isn't that uncommon, by the way, it would essentially be a sort of arbitrary-precision numeric type.

The idea is not to go on forever retaining all the accuracy after every operation. It's just that my original types have limited accuracy, and I want to control when to lose it. In practice the operations would return the highest accuracy, and I have to control when to throw away bits. While this may not be practical when prototyping, it might be useful to code a processing pipeline that makes extensive use of 8-bit and 16-bit SIMD operations, for instance. My N1f7 example is from a real application, and right now if I want to code an optimized fixed-precision pipeline I'll have to make the low-level computations by myself.

I'm curious to experimenting with this accuracy-preserving type, although I'm not sure it's really interesting to many people. Is this something you imagine could ever be a part of this package, or is it better if I develop something separate?

[kimikage]

I cannot determine if it could be part of this package without a more concrete example.

There are two things to note. Most users only use N0f8, some use N0f16, and the other types are rarely used. Also, 2^f - 1 is very awkward when considering the accuracy (not precision).

[timholy]

This isn't that uncommon, by the way, it would essentially be a sort of arbitrary-precision numeric type.

which are painfully slow:

julia> @btime x*y setup=(x = rand(Int); y = rand(Int))
  1.239 ns (0 allocations: 0 bytes)
-5710601054559440257

julia> @btime x*y setup=(x = big(rand(Int)); y = big(rand(Int)))
  59.760 ns (3 allocations: 48 bytes)
-28182581540266793182579202048252613136

We need to use hardware-supported arithmetic wherever possible.

I have to control when to throw away bits

Yep. If you need a workaround now...often you can do that by manually converting one value at the outset to the final precision you want to achieve. That's why I recommended float, but you could also do this:

julia> x = rand(N0f8)
0.514N0f8

julia> w = rand(N1f7)
0.339N1f7

julia> N1f15(x)*w
0.17392N1f15

and the other types are rarely used

I agree with that overall. In the past my lab used N2f14 a lot because we had a 14-bit scientific camera, but all of our cameras in current use are 16-bit. That said, there are some attractive cameras that are 10-bit, so some of those other types could come back into usage.

[kimikage]

It would be good to have some specific and practical model cases or benchmark tests.
In my use cases, I rarely run into "unintended" overflow problems. I may be in a too low place to realize that I am falling into a pit. 😅

We need to consider a little more about the workflow of the development of the program than the code itself. I think the workflow will also change depending on the nature of the images and the purpose (which is more important, the results or the methods?).

[johnnychen94]

reduce(+, A)/length(A) , mapreduce(abs2, +, A) and almost all the distances in Distances.jl if you try to implement it from scratch 😂

[kimikage]

Is this really a practical case for the first two? I think it's clear that they don't work well "as is" in many cases, even if we automatically widen the 8-bit type to the 16-bit type. (Of course, as you know, there is a difference between reduce and sum. https://docs.julialang.org/en/v1/base/collections/#Base.sum)

I am not saying that they are trivial. I believe that the solutions will depend on whether they are practical or not. One of the reasons I mentioned workflows is to separate cases of coding errors (requiring code fixes outside of FixedPointNumbers) from the rest. Since Julia's static analysis tools are still underdeveloped, an effective countermeasure against coding errors is the checked arithmetic. If the one-time widening arithmetic is the default, the cases of coding errors are more aggravated. 😱 This is why I don't like the "implicit" widening arithmetic.

[timholy]

julia> using TestImages, ImageView, ImageCore

julia> img = testimage("cameraman");

julia> scaling = LinRange(0.9, 1.1, size(img, 1));   # let's simulate a change in lighting

julia> imgs = Gray{N0f8}.(clamp.(scaling .* img, 0, 1));   # in the real world this would be a second image, so it would also have Gray{N0f8} eltype

julia> imshow(img - imgs)
Dict{String, Any} with 4 entries:
  "gui"         => Dict{String, Any}("window"=>GtkWindowLeaf(name="", parent, width-request=-1, height-request=-1, visible=TRUE, sensitive=TRUE, app-paintable=FALSE, can-focus=FALSE, has-focus=FALSE, is-…
  "roi"         => Dict{String, Any}("redraw"=>122: "map(clim-mapped image, input-32)" = nothing Nothing , "zoomregion"=>86: "input-31" = ZoomRegion{RInt64}(XY(1..512, 1..512), XY(1..512, 1..512)) ZoomRe…
  "annotations" => 88: "input-32" = Dict{UInt64, Any}() Dict{UInt64, Any} 
  "clim"        => 87: "CLim" = CLim{N0f8}(0.0, 1.0) CLim{N0f8} 

julia> imshow(float.(img) .- float.(imgs))
Dict{String, Any} with 4 entries:
  "gui"         => Dict{String, Any}("window"=>GtkWindowLeaf(name="", parent, width-request=-1, height-request=-1, visible=TRUE, sensitive=TRUE, app-paintable=FALSE, can-focus=FALSE, has-focus=FALSE, is-…
  "roi"         => Dict{String, Any}("redraw"=>130: "map(clim-mapped image, input-44)" = nothing Nothing , "zoomregion"=>94: "input-43" = ZoomRegion{RInt64}(XY(1..512, 1..512), XY(1..512, 1..512)) ZoomRe…
  "annotations" => 96: "input-44" = Dict{UInt64, Any}() Dict{UInt64, Any} 
  "clim"        => 95: "CLim" = CLim{Float32}(-0.0901961, 0.0745098) CLim{Float32} 

First image:

Second image:

The first is completely useless.

This is a day-in, day-out operation in scientific image processing, where you often want to know what changed over time or how your image registration algorithm changed the alignment of two images.

[kimikage]

I'm skeptical that overflow still can be a practical issue, even though that is a day-in, day-out operation. 😕
My daily task is as follows, though. 😭

(@v1.5) pkg> git checkout master
ERROR: Could not determine command

I think it's a good benchmark, but I think the latter is scientifically "bogus". We should take account of the extra color mapping. It must have more costs than the subtraction. (cf. #143 (comment))

[timholy]

I don't understand your meaning in your previous post, can you clarify? (Neither the git command example nor the extra color mapping. Sorry!)

[kimikage]

Sorry for the confusing joke. I was trying to say that people make mistakes on a daily basis, even if they know it's a problem, and that they rarely take steps to prevent it from happening again. And there is no doubt that the fool-proof "checking" is effective and annoying. 😄

I'm highlighting the point that imshow (and its back-end) detects the maximum and minimum pixel values, determines the color scale accordingly, and finally generates the "RGB" image in VRAM. If this is a common case, then the conclusion would be that introducing checked_sub would have little impact on performance. 😅

[johnnychen94]

actually I don’t think it’s a mistake, because it’s working well in MATLAB (not exactly, but kind of) :P How many such “bugs” do you expect people meet until they seriously check this behavior in their code? I would say 0, at most 1. But should never be a daily issue. We’re discussing this behavior at a state of knowing this and being capable of debugging it instantly. But there are many non-developer people out there would be bitten by this and find it difficult to realize the actual issue. There will be more when Images get stable and popular. This is mostly of my personal experience and can be biased. Johnny Chen 陈久宁 School of Mathematical Sciences, East China Normal University, Shanghai 华东师范大学数学科学学院 2020年9月16日 +0800 19:33 kimikage <[email protected]>，写道：

…

Sorry for the confusing joke. I was trying to say that people make mistakes on a daily basis, even if they know it's a problem, and that they rarely take steps to prevent it from happening again. And there is no doubt that the fool-proof "checking" is effective and annoying. 😄 I'm highlighting the point that imshow (and its back-end) detects the maximum and minimum pixel values, determines the color scale accordingly, and finally generates the "RGB" image in VRAM. If this is a common case, then the conclusion would be that introducing checked_sub would have little impact on performance. 😅 — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[kimikage]

I'm not looking for clarity on how often overflows happen and how many people they annoy, but I'm looking for clarity on where in the workflow they occur. I don't want to use JuMP to propose a solution that maximizes overall satisfaction, but rather to find the best practices for each case and build tools and documentation to support each practice.

[johnnychen94]

Speaking of the best practice, I think Tim's example has shown it quite clearly, that the safest strategy right now is to promote to a computation type; float in our case. Tim's proposal goes one more step in this proposal, to make this "promotion" happen automatically whenever necessary with the least cost.

Let's take a non-concrete example. In practice, there might be codes that we can't control of, say, function f(arr) developed by another community that never considered to support FixedPointNumbers. It is understandable for them to not promote the element types. How could we make use of this f without worrying about the overflow issues? The best and perhaps the only strategy is to add a wrapper to it:

f(arr::AbstractArray{<:Normed}) = f(MappedArrays.of_eltype(float(eltype(arr)), arr))

Could we avoid this manually written wrapper? IIUC Tim's proposal makes it possible.

[kimikage]

That is off-topic, i.e. more "short-sighted" than what we're discussing here, but I posted a topic on the Discourse about the default arithmetic.
https://discourse.julialang.org/t/rfc-what-should-the-arithmetic-within-the-fixedpointnumbers-be/46697
I'm not a great communicator of the written word, so I would appreciate it if you could help me get more people interested in the overflow issue.

[kimikage]

Tim's proposal goes one more step in this proposal, to make this "promotion" happen automatically whenever necessary with the least cost.

Ideally, yes, but in reality, this is not the case. The information we can get from the types (and operations) is very limited. (The array, or image, is a bit more informative about the overflow safety in the folding operation, though).

[nlw0]

I know my example is maybe too niche. I think it would be good to talk about specific applications, though, so we can better understand what to expect from this library, and whether there's a potential to create a new library perhaps.

Also let me make it clear I understand most users only want this library for storing and reading data, and talking about how to operate is already "too niche" in a way.

One idea of a prototypical application I have right now is linear filters. Of course you should just move to floating point if you can. There are situations filtering audio or image signals where fixed point can be good enough, while delivering significant speed gains. Has anyone ever tried to implement a filter using this library? Are there any shortcomings? I would love to collaborate with more people trying that out and seeing how it goes, maybe even actually trying to interface with the DSP and SIMD libraries. Any volunteers?

[kimikage]

This is not a decision, but I have planned to support the 3-arg methods which can specify the destination type manually. This plan leaves the user in control of overflow, accuracy and precision for now.

I've used FIRs and IIRs with few taps, but the filter design itself was done in floating point numbers, so I didn't find it particularly inconvenient. Of course, that's just the case of my application, and I agree that there is room for improvement there.

My concern is the rounding mode. A variety of rounding modes are used in fixed-point arithmetic (although RoundNearest is rather minor), but their support is daunting.

I still don't fully understand what you want to do, so I think it would be good to flesh it out in another issue.

Edit:
In addition, I guess @timholy was planning to adapt some of the JuliaImages ecosystem to LoopVectorization. I don't know the technical details of that, but it would be helpful in terms of SIMD.