code

REQUIRE

@@ -1 +1,2 @@
 julia 0.6
+StaticArrays

src/DynamicApproximations.jl

@@ -1,5 +1,10 @@
 module DynamicApproximations
 
+using StaticArrays
+
+include("types/LList.jl")
+include("types/QuadraticPolynomial.jl")
+include("types/PiecewiseQuadratic.jl")
 # package code goes here
 
 end # module

src/types/LList.jl

@@ -0,0 +1,58 @@
+type LList{T}
+    val::T
+    next::LList{T}
+    LList{T}() where T = new() #Hack for empty list
+    LList{T}(val::T) where T = (l = new(val); l.next = l)
+    LList{T}(val::T, l::LList{T}) where T = new(val, l)
+end
+LList{T}(val::T) = LList{T}(val)
+LList{T}(val::T, l::LList{T}) = LList{T}(val, l)
+
+llist{T}(val::T) = LList(val)
+llist{T}(val::T, args...) = LList(val, llist(args...))
+#Hack for empty set
+Base.isempty(l::LList) = !isdefined(l, :next)
+
+function Base.show{T}(io::IO, l::LList{T})
+    print(io, "llist(")
+    if !isempty(l)
+        show(io, l.val)
+        if !islast(l)
+            for val in l.next
+                print(io, ",")
+                show(io, val)
+            end
+        end
+    end
+    print(io, ")")
+end
+
+islast(l::LList) = l === l.next || isempty(l.next)
+
+"""
+Insert `val` in `l=llist(v1,v2,...,vn)` so that
+`l=llist(v1,val,v2,...,vn)`.
+If `l` is empty then `l` becomes `llist(val)`.
+"""
+function insertafter!{T}(val::T, l::LList{T})
+    if isempty(l)
+        l.val = val
+        l.next = l
+    else
+        l2 = islast(l) ? LList{T}(val) : LList{T}(val, l.next)
+        l.next = l2
+    end
+    return l
+end
+
+Base.start{T}(l::LList{T}) = (l, isempty(l))
+Base.done{T}(l::LList{T}, state::Tuple{LList{T},Bool}) = state[2]
+Base.next{T}(l::LList{T}, state::Tuple{LList{T},Bool}) = (state[1].val, (state[1].next, islast(state[1])))
+
+function Base.length(l::LList)
+    n = 0
+    for i in l
+        n += 1
+    end
+    n
+end

src/types/PiecewiseQuadratic.jl

@@ -0,0 +1,28 @@
+struct QuadraticOnInterval{T}
+    p::QuadraticPolynomial{T}
+    endpoint::T
+end
+
+struct PiecewiseQuadratic{T}
+    qlist::LList{QuadraticOnInterval{T}}
+end
+
+PiecewiseQuadratic{T}(p::QuadraticPolynomial{T}) = PiecewiseQuadratic{T}(llist(QuadraticOnInterval(p,typemax(T))))
+
+"""
+    min!(q1::PiecewiseQuadratic{T}, q2::QuadraticOnInterval{T})
+Update `q1(x)` on the interval `I` in `q2` to be `q1(x) := min(q1(x),q2(x)) ∀ x∈I`
+"""
+function min!{T}(q1::PiecewiseQuadratic{T}, q2::QuadraticOnInterval{T})
+    #TODO
+end
+
+
+Λs = Array{PiecewiseQuadratic,1}(len1,len2)
+N = 10
+M = 4
+for m = M:-1:1
+    for i = (N-m):-1:1
+        Λim = PiecewiseQuadratic(QuadraticPolynomial(Inf,Inf,Inf))
+        for ip = (i+1):(N-m+1)
+            update!(Λim, ts[ip], Λ[m,ip], l[i,ip])

src/types/QuadraticPolynomial.jl

@@ -0,0 +1,51 @@
+struct QuadraticPolynomial{T<:Real}
+    a::T
+    b::T
+    c::T
+end
+
+function QuadraticPolynomial{T}(a::T,b::T,c::T)
+    @assert a ≥ 0
+    QuadraticPolynomial{T}(a,b,c)
+end
+
+"""
+    n, intersec = intersections{T}(p1::QuadraticPolynomial{T},p2::QuadraticPolynomial{T})
+    returns number of intersections `n` and a 2-vector `intersec` containing the intersection points
+    in the first `n` elements.
+"""
+function intersections{T}(p1::QuadraticPolynomial{T},p2::QuadraticPolynomial{T})
+    a = p1.a-p2.a
+    b = p1.b-p2.b
+    c = p1.c-p2.c
+    n = 0
+    intersec = [T(NaN), T(NaN)]
+    if a == 0
+        #Numerical problems if division, 0 or 1 intersections
+        if b == 0
+            # 0 intersections
+            #return 0, [T(NaN), T(NaN)]
+        else
+            # 1 intersection
+            n = 1
+            intersec[1] = -c/b
+            #return 1, [-c/b, T(NaN)]
+        end
+    else
+        # 0 or 2 intersections, possibly in same point
+        r = b^2 - 4c*a
+        if r < 0
+            # 0 intersections
+            #return 0, [T(NaN), T(NaN)]
+        else
+            # 2 intersections
+            r = sqrt(r)/(2a)                 #sqrt(b^2-4ca)/2a
+            v = -b/(2a)
+            n = 2
+            intersec[1] = v-r
+            intersec[2] = v+r
+            #return 2, [v-r, v+r]
+        end
+    end
+    return n, intersec
+end