Julia Season of Contributions Blog#1

My project titled 'Standard Compliant IntervalArithmetic Library' got selected in Julia Season of Contribution under the mentorship of David P. Sanders. My project is under JuliaIntervals to develop IntervalArithmetic.jl and make it compliant with the IEEE1788-2015, a paper which provides standards for interval computations. Basically my project comprised of including the missing functions and improving the existing functions upto the standards.

Community Bonding

I started with reading the standard for the implementation of the Reduction and String Parsing functions. In the meantime, I also started to solve some issues, namely cospi/sinpi stack overflow and broken rational power functions. I learned how to write regular expressions for the string parsing functions.

Starting with the project

During the end of the community boding period I started with the implementation of the functions planned in my proposal. I had reduction functions and the string parsing functions in my first two weeks. Lets see how did they go in the implementation section.

Implementation

Sinpi and Cospi functions

sinPi and cosPi are defined as sin(πx) and cos(πx) respectively with domain as R and range equal to [-1, 1]. These functions should return a range containing the sine and cosine of values in the interval πx. I used the implementation of sin and cos functions to implement sinPi and cosPi functions in IntervalArithmetic.jl.

import Pkg
Pkg.add("IntervalArithmetic")
using IntervalArithmetic
function sinpi(a::Interval{T}) where T
    isempty(a) && return a
    w = a * pi_interval(T)
    return sin(w)
end
sinpi(0.5 .. 0.75)

2.1s

Julia

[0.707106, 1]

Fixing Rational Power

The function ^(a :: Interval, x :: Rational) computes the range of values in the interval raised to the power x. This issue was not straight as it was ambiguous on some platform. First I started with using the power functions by type casting x into BigFloat. The added tests worked well but it failed for some existing tests. Then I moved onto compute the qth root raised to the power p where p is numerator and q is the denominator of x. I calculated the lower and upper bounds for qth root but I could not get the correct rounding as per the standard. Then I used ccall function to wrap nthroot function in GNU MPFR library to ensure correct rounded results.

function ^(a::Interval{T}, x::Rational) where T
  	if x > 0
        low = a.lo ^ BigFloat(x)
        high = a.hi ^ BigFloat(x)
        if isinteger(high) && isinteger(low)
            return Interval(low, high)
        end
        return ^(a, BigFloat(x))
    end
    if x < 0
        return inv(^(a, -x))
    end
end

0.6s

Julia

^ (generic function with 1 method)

function ^(a::Interval{T}, x::Rational) where T
   p = x.num
   q = x.den
   if x > 0
        low = a.lo ^ BigFloat(1//q)
        high = a.hi ^ BigFloat(1//q)
        isinteger(high) && isinteger(low) && return ^(Interval(low, high) , p)
        b = @interval(low, high)
        return b ^ p
   end

Shift+Enter to run

Julia

0.6s

Julia

import Base.MPFR: MPFRRoundingMode
import Base.MPFR: MPFRRoundUp, MPFRRoundDown, MPFRRoundNearest, MPFRRoundToZero, MPFRRoundFromZero
function ^(a::Interval{T}, x::Rational) where T
	  domain = Interval{T}(0, Inf)

    p = x.num
    q = x.den

    isempty(a) && return emptyinterval(a)
    x == 0 && return one(a)
		if a == zero(a)
    x > zero(x) && return zero(a)
        return emptyinterval(a)
    end
  	x == (1//2) && return sqrt(a)

    if x >= 0
        if a.lo ≥ 0
            isinteger(x) && return a ^ Int64(x)
            a = @biginterval(a)
            ui = convert(Culong, q)
            low = BigFloat()
            high = BigFloat()
            ccall((:mpfr_rootn_ui , :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , a.lo , ui, MPFRRoundDown)
            ccall((:mpfr_rootn_ui , :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , a.hi , ui, MPFRRoundUp)
            b = @interval(low , high)
            return b^p
        end
    	if a.lo < 0 && a.hi ≥ 0
            isinteger(x) && return a ^ Int64(x)
            a = a ∩ Interval{T}(0, Inf)
            a = @biginterval(a)
            ui = convert(Culong, q)
            low = BigFloat()
            high = BigFloat()
            ccall((:mpfr_rootn_ui , :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , a.lo , ui, MPFRRoundDown)
            ccall((:mpfr_rootn_ui , :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , a.hi , ui, MPFRRoundUp)
            b = @interval(low , high)
            return b^p
        end

        if a.hi < 0
            isinteger(x) && return a ^ Int64(x)
            return emptyinterval(a)
        end
    end

    if x < 0
        return inv(a^(-x))
    end
end

0.6s

Julia

^ (generic function with 1 method)

Reduction functions

We have four operations on vectors, namely sum, dot, sumAbs and sumSquare. This implementation is currently under review. There are issues related to the precision of BigFloat and generating different tests for different rounding modes. The standard states that these functions should return the correctly rounded results hence I used BigFloat and then rounded it to return Float64 value.

function vector_sum(v :: Vector{T} , r :: RoundingMode) where T
  sum1 = sum(BigFloat.(v))
  return (sum1, r)
end

function vector_dot(v :: Vector{T}, u :: Vector{T}, r :: RoundingMode) where T
  sum1 = sum(BigFloat.(v) .* BigFloat.(u))
  return (sum1, r)
end

function vector_sumAbs(v :: Vector{T}, r :: RoundingMode) where T
  sum1 = sum(BigFloat.(abs.(v)))
  return (sum1, r)
end

function vector_sumSquare(v :: Vector{T}, r :: RoundingMode) where T
  return vector_dot(v, v, r)
end

0.5s

Julia

(15.0, RoundingMode{:Up}())

String Parsing

I had to implement mostly the uncertain form(like 2.4?1) conversion. Moreover it required recognition of regular expressions for different types of strings.

0.6s

Julia

function parse(::Type{Interval{T}}, s::AbstractString) where T
	m = match(r"(\-?\d*\.?\d*)\?\?([ud]?)", s)   # match with string of form "34??"

            if m!= nothing
                if m.captures[2] == ""    # strings of the form "10??"
                    return entireinterval(T)
                end
                if m.captures[2] == "u"   # strings of the form "10??u"
                    lo = parse(Float64, m.captures[1])
                    hi = Inf
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end
                if m.captures[2] == "d"   # strings of the form "10??d"
                    lo = -Inf
                    hi = parse(Float64, m.captures[1])
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end
            end

            m = match(r"(\-?\d*)\.?(\d*)\?(\d*)([ud]?)(e-?\d*)?", s) # match with strings like "3.56?1u" or "3.56?1e2"

            if m != nothing  # matched
                if m.captures[3] != "" && m.captures[4] == "" && m.captures[5] == nothing # string of form "3.4?1" or "10?2"
                    d = length(m.captures[2])
                    x = parse(Float64, m.captures[3] * "e-$d")
                    n = parse(Float64, m.captures[1]*"."*m.captures[2])
                    lo = n - x
                    hi = n + x
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end

                if m.captures[3] == "" && m.captures[4] == "" && m.captures[5] == nothing  # strings of the form "3.46?"
                    d = length(m.captures[2])
                    x = parse(Float64, "0.5" * "e-$d")
                    n = parse(Float64, m.captures[1]*"." * m.captures[2])
                    lo = n - x
                    hi = n + x
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end

                if m.captures[3] == "" && m.captures[4] == "" && m.captures[5] != nothing  # strings of the form "3.46?e2"
                    d = length(m.captures[2])
                    x = parse(Float64, "0.5" * "e-$d")
                    n = parse(Float64, m.captures[1]*"." * m.captures[2])
                    lo = parse(Float64, string(n-x) * m.captures[5])
                    hi = parse(Float64, string(n+x) * m.captures[5])
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end

                if m.captures[3] == "" && m.captures[4] == "u" && m.captures[5] == nothing # strings of the form "3.4?u"
                    d = length(m.captures[2])
                    x = parse(Float64, "0.5" * "e-$d")
                    n = parse(Float64, m.captures[1]*"."*m.captures[2])
                    lo = n
                    hi = n+x
                    interval = eval(make_interval(T, lo, [hi]))
                    return interval
                end

                if m.captures[3] == "" && m.captures[4] == "d" && m.captures[5] == nothing # strings of the form "3.4?d"
                        d = length(m.captures[2])
                        x = parse(Float64, "0.5" * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = n - x
                        hi = n
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                end

                if m.captures[3] != ""
                    if m.captures[4] == "u" && m.captures[5] != nothing    #strings of the form "3.46?1u"
                        d = length(m.captures[2])
                        x = parse(Float64, m.captures[3] * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = parse(Float64, string(n) * m.captures[5])
                        hi = parse(Float64, string(n+x) * m.captures[5])
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                    end

                    if m.captures[4] == "u" && m.captures[5] == nothing    #strings of the form "3.46?1u"
                        d = length(m.captures[2])
                        x = parse(Float64, m.captures[3] * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = n
                        hi = n+x
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                    end

                    if m.captures[4] == "d" && m.captures[5] != nothing  #strings of the form "3.46?1d"
                        d = length(m.captures[2])
                        x = parse(Float64, m.captures[3] * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = parse(Float64, string(n-x) * m.captures[5])
                        hi = parse(Float64, string(n) * m.captures[5])
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                    end

                    if m.captures[4] == "d" && m.captures[5] == nothing  #strings of the form "3.46?1d"
                        d = length(m.captures[2])
                        x = parse(Float64, m.captures[3] * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = n-x
                        hi = n
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                    end

                    if m.captures[4] == "" && m.captures[5] != nothing    # strings of the form "3.56?1e2"
                        d = length(m.captures[2])
                        x = parse(Float64, m.captures[3] * "e-$d")
                        n = parse(Float64, m.captures[1]*"."*m.captures[2])
                        lo = parse(Float64, string(n-x)*m.captures[5])
                        hi = parse(Float64, string(n+x)*m.captures[5])
                        interval = eval(make_interval(T, lo, [hi]))
                        return interval
                    end
                end

            m = match(r"(-?\d*\.?\d*)", s)  # match strings of form "1" and "2.4"

            if m != nothing
                lo = parse(Float64, m.captures[1])
                hi = lo
            end

            if m == nothing
                throw(ArgumentError("Unable to process string $s as interval"))
            end

            interval = eval(make_interval(T, lo, [hi]))
            return interval
        end
end

0.6s

Julia

parse (generic function with 1 method)

What next?

I am following my project timeline and will try to do so during whole project. I have IntervalToText conversion next on my list and some other functions. I will also try to fix some issues apart from these. I must say its been great experience with the Julia community till now and I am looking forward to the completion of my project.