Rounding Algorithms from Property Based Testing

Steffen Forkmann recently posted a blog about incorrect rounding in a twitter poll, and how coding a rounding strategy is a hard problem. This got me thinking about how many correct rounding algorithms there are.

In these type of problems it is important to look at what properties the algorithm should have. Property based testing is a great tool when doing this.

The key property required for a fair rounding algorithm is that rounded values increase with the weights. It doesn't make sense for a lower weight to have a greater rounded value. Symmetry in the results for negative weights and negative value to be distributed are also important. This can easily be achieved by mapping from the positive results, but a robust algorithm shouldn't need to do this.

In the blog it was proposed that adjusting the largest weight would work, but in general this can only work when the rounding needs a positive adjustment due to the increasing with weight property. For negative adjustments the smallest non-zero weight would need to be adjusted. It may also be unfair to adjust these values if they already have a large rounding error.

Error minimisation algorithm

This algorithm minimises the absolute and relative rounding errors. I normally apply absolute then relative but the reverse order also works and may be more correct for certain problems.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
/// Distribute integer n over an array of weights
let distribute n weights =
    let wt = Array.sum weights
    if Array.isEmpty weights
        || Seq.maxBy abs weights * float(abs n+1) |> abs
           >= abs wt * float Int32.MaxValue then None
    else
        let inline round f = int(if f<0.0 then f-0.5 else f+0.5)
        let ns = Array.map ((*)(float n / wt) >> round) weights

        let inline absoluteError wi ni d =
            let wc = 0.5 * float d * wt / float n
            let wni = float ni * wt / float n
            if wc > 0.0 then min (wni-wi) 0.0 + wc
            else min (wi-wni) 0.0 - wc
            
        let d = n - Array.sum ns
        for __ = 1 to abs d do
            let _,_,_,i =
                Seq.mapi2 (fun i wi ni ->
                    let absErr = absoluteError wi ni (sign d)
                    let relErr = absErr / abs wi
                    let weight = if wt > 0.0 then -wi else wi
                    absErr, relErr, weight, i
                ) weights ns |> Seq.min
            ns.[i] <- ns.[i] + sign d
        Some ns

Testing

The twitter tricky test below is interesting. It is not clear which values should be adjusted down. Neither the largest or smallest weights look like good candidates. The error minimisation algorithm sensibly selects the second largest weight and keeps the correct order.

Testing also shows the algorithm is sensitive to weights that are close together and issues with machine precision. This directs how the error function needs to calculate.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
42: 
43: 
44: 
45: 
46: 
47: 
48: 
49: 
50: 
51: 
52: 
53: 
54: 
55: 
56: 
57: 
58: 
59: 
60: 
61: 
62: 
63: 
64: 
65: 
66: 
67: 
68: 
69: 
70: 
71: 
72: 
73: 
74: 
75: 
76: 
77: 
78: 
79: 
80: 
81: 
82: 
83: 
84: 
85: 
86: 
87: 
88: 
89: 
90: 
91: 
92: 
93: 
94: 
95: 
96: 
let roundingTests =
    testList "rounding" [

        test "empty" {
            let r = distribute 1 [||]
            Expect.equal r None "empty"
        }

        test "n zero" {
            let r = distribute 0 [|406.0;348.0;246.0;0.0|]
            Expect.equal r (Some [|0;0;0;0|]) "zero"
        }

        test "twitter" {
            let r = distribute 100 [|406.0;348.0;246.0;0.0|]
            Expect.equal r (Some [|40;35;25;0|]) "40 etc"
        }

        test "twitter n negative" {
            let r = distribute -100 [|406.0;348.0;246.0;0.0|]
            Expect.equal r (Some [|-40;-35;-25;0|]) "-40 etc"
        }

        test "twitter weights negative" {
            let r = distribute 100 [|-406.0;-348.0;-246.0;-0.0|]
            Expect.equal r (Some [|40;35;25;0|]) "40 etc"
        }

        test "twitter both negative" {
            let r = distribute -100 [|-406.0;-348.0;-246.0;-0.0|]
            Expect.equal r (Some [|-40;-35;-25;0|]) "-40 etc"
        }

        test "twitter tricky" {
            let r = distribute 100 [|404.0;397.0;57.0;57.0;57.0;28.0|]
            Expect.equal r (Some [|40;39;6;6;6;3|]) "o no"
        }

        test "negative example" {
            let r1 = distribute 42 [|1.5;1.0;39.5;-1.0;1.0|]
            Expect.equal r1 (Some [|2;1;39;-1;1|]) "2 etc"
            let r2 = distribute -42 [|1.5;1.0;39.5;-1.0;1.0|]
            Expect.equal r2 (Some [|-2;-1;-39;1;-1|]) "-2 etc"
        }

        testProp "ni total correctly" (fun n (Gen.RationalFloats ws) ->
            distribute n ws
            |> Option.iter (fun ns -> Expect.equal (Array.sum ns) n "sum ns = n")
        )

        testProp "negative n returns opposite of positive n" (
            fun n (Gen.RationalFloats ws) ->
                let r1 = distribute -n ws |> Option.map (Array.map (~-))
                let r2 = distribute n ws
                Expect.equal r1 r2 "r1 = r2"
        )

        testProp "increase with weight" (fun n (Gen.RationalFloats ws) ->
            let d = if Seq.sum ws > 0.0 <> (n>0) then -1 else 1
            distribute n ws
            |> Option.iter (
                   Seq.map ((*)d)
                >> Seq.zip ws
                >> Seq.sort
                >> Seq.pairwise
                >> Seq.iter (fun (ni1,ni2) ->
                    Expect.isLessThanOrEqual ni1 ni2 "ni1 <= ni2")
            )
        )

        testProp "smallest error" (fun n (Gen.RationalFloats ws) randomChanges ->
            distribute n ws
            |> Option.iter (fun ns ->
                let totalError ns =
                    let wt = Seq.sum ws
                    Seq.map2 (fun ni wi ->
                        float ni * wt / float n - wi |> abs
                    ) ns ws
                    |> Seq.sum

                let errorBefore = totalError ns

                let l = Array.length ns
                List.iter (fun (i,j) ->
                    ns.[abs i % l] <- ns.[abs i % l] - 1
                    ns.[abs j % l] <- ns.[abs j % l] + 1
                ) randomChanges

                let errorAfter = totalError ns

                Expect.floatLessThanOrClose
                    Accuracy.veryHigh errorBefore errorAfter "before <= after"
                )
        )

    ]

Conclusion

This is an example of how property based testing can actually help in algorithm design. It gives example failing cases that can steer you to a better solution.

I have seen this problem in order management systems where orders for a number of shares are to be allocated across a number of portfolios. The buy and sell orders have a number of partial fills, but ultimately everything needs to add up in a consistent and robust way.

Error minimisation is the best rounding algorithm I have found. By construction it is also the one with the smallest total rounding error.

I don't think there is any other general purpose algorithm that can satisfy these properties as well.

module Rounding
namespace System
val distribute : n:int -> weights:float [] -> int [] option

Full name: Rounding.distribute


 Distribute integer n over an array of weights
val n : int
val weights : float []
val wt : float
type Array =
  member Clone : unit -> obj
  member CopyTo : array:Array * index:int -> unit + 1 overload
  member GetEnumerator : unit -> IEnumerator
  member GetLength : dimension:int -> int
  member GetLongLength : dimension:int -> int64
  member GetLowerBound : dimension:int -> int
  member GetUpperBound : dimension:int -> int
  member GetValue : [<ParamArray>] indices:int[] -> obj + 7 overloads
  member Initialize : unit -> unit
  member IsFixedSize : bool
  ...

Full name: System.Array
val sum : array:'T [] -> 'T (requires member ( + ) and member get_Zero)

Full name: Microsoft.FSharp.Collections.Array.sum
val isEmpty : array:'T [] -> bool

Full name: Microsoft.FSharp.Collections.Array.isEmpty
module Seq

from Microsoft.FSharp.Collections
val maxBy : projection:('T -> 'U) -> source:seq<'T> -> 'T (requires comparison)

Full name: Microsoft.FSharp.Collections.Seq.maxBy
val abs : value:'T -> 'T (requires member Abs)

Full name: Microsoft.FSharp.Core.Operators.abs
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>
type Int32 =
  struct
    member CompareTo : value:obj -> int + 1 overload
    member Equals : obj:obj -> bool + 1 overload
    member GetHashCode : unit -> int
    member GetTypeCode : unit -> TypeCode
    member ToString : unit -> string + 3 overloads
    static val MaxValue : int
    static val MinValue : int
    static member Parse : s:string -> int + 3 overloads
    static member TryParse : s:string * result:int -> bool + 1 overload
  end

Full name: System.Int32
field int.MaxValue = 2147483647
union case Option.None: Option<'T>
val round : (float -> int)
val f : float
Multiple items
val int : value:'T -> int (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.int

--------------------
type int = int32

Full name: Microsoft.FSharp.Core.int

--------------------
type int<'Measure> = int

Full name: Microsoft.FSharp.Core.int<_>
val ns : int []
val map : mapping:('T -> 'U) -> array:'T [] -> 'U []

Full name: Microsoft.FSharp.Collections.Array.map
val absoluteError : (float -> 'a -> 'b -> float) (requires member op_Explicit and member op_Explicit)
val wi : float
val ni : 'a (requires member op_Explicit)
val d : 'b (requires member op_Explicit)
val wc : float
val wni : float
val min : e1:'T -> e2:'T -> 'T (requires comparison)

Full name: Microsoft.FSharp.Core.Operators.min
val d : int
val __ : int
val i : int
val mapi2 : mapping:(int -> 'T1 -> 'T2 -> 'U) -> source1:seq<'T1> -> source2:seq<'T2> -> seq<'U>

Full name: Microsoft.FSharp.Collections.Seq.mapi2
val ni : int
val absErr : float
val sign : value:'T -> int (requires member get_Sign)

Full name: Microsoft.FSharp.Core.Operators.sign
val relErr : float
val weight : float
val min : source:seq<'T> -> 'T (requires comparison)

Full name: Microsoft.FSharp.Collections.Seq.min
union case Option.Some: Value: 'T -> Option<'T>
Multiple items
union case RationalFloats.RationalFloats: float [] -> RationalFloats

--------------------
type RationalFloats = | RationalFloats of float []

Full name: Rounding.Gen.RationalFloats
val rationalFloats : obj

Full name: Rounding.Gen.rationalFloats
val fraction : (int -> int -> int -> float)
val a : int
val b : int
val c : int
val private config : obj

Full name: Rounding.config
val typeof<'T> : Type

Full name: Microsoft.FSharp.Core.Operators.typeof
module Gen

from Rounding
Multiple items
union case Gen.RationalFloats.RationalFloats: float [] -> Gen.RationalFloats

--------------------
type RationalFloats = | RationalFloats of float []

Full name: Rounding.Gen.RationalFloats
val testProp : name:'a -> 'b

Full name: Rounding.testProp
val name : 'a
val ptestProp : name:'a -> 'b

Full name: Rounding.ptestProp
val ftestProp : stdgen:'a -> name:'b -> 'c

Full name: Rounding.ftestProp
val stdgen : 'a
val name : 'b
val roundingTests : obj

Full name: Rounding.roundingTests
module Option

from Microsoft.FSharp.Core
val iter : action:('T -> unit) -> option:'T option -> unit

Full name: Microsoft.FSharp.Core.Option.iter
val map : mapping:('T -> 'U) -> option:'T option -> 'U option

Full name: Microsoft.FSharp.Core.Option.map
val sum : source:seq<'T> -> 'T (requires member ( + ) and member get_Zero)

Full name: Microsoft.FSharp.Collections.Seq.sum
val map : mapping:('T -> 'U) -> source:seq<'T> -> seq<'U>

Full name: Microsoft.FSharp.Collections.Seq.map
val zip : source1:seq<'T1> -> source2:seq<'T2> -> seq<'T1 * 'T2>

Full name: Microsoft.FSharp.Collections.Seq.zip
val sort : source:seq<'T> -> seq<'T> (requires comparison)

Full name: Microsoft.FSharp.Collections.Seq.sort
val pairwise : source:seq<'T> -> seq<'T * 'T>

Full name: Microsoft.FSharp.Collections.Seq.pairwise
val iter : action:('T -> unit) -> source:seq<'T> -> unit

Full name: Microsoft.FSharp.Collections.Seq.iter
val map2 : mapping:('T1 -> 'T2 -> 'U) -> source1:seq<'T1> -> source2:seq<'T2> -> seq<'U>

Full name: Microsoft.FSharp.Collections.Seq.map2
val length : array:'T [] -> int

Full name: Microsoft.FSharp.Collections.Array.length
Multiple items
module List

from Microsoft.FSharp.Collections

--------------------
type List<'T> =
  | ( [] )
  | ( :: ) of Head: 'T * Tail: 'T list
  interface IEnumerable
  interface IEnumerable<'T>
  member GetSlice : startIndex:int option * endIndex:int option -> 'T list
  member Head : 'T
  member IsEmpty : bool
  member Item : index:int -> 'T with get
  member Length : int
  member Tail : 'T list
  static member Cons : head:'T * tail:'T list -> 'T list
  static member Empty : 'T list

Full name: Microsoft.FSharp.Collections.List<_>
val iter : action:('T -> unit) -> list:'T list -> unit

Full name: Microsoft.FSharp.Collections.List.iter