Merge branch 'master' of github.com:JuliaLang/julia

extras/ode.jl

@@ -51,7 +51,7 @@ function ode23(F::Function, tspan::AbstractVector, y_0::AbstractVector)
     # Compute initial step size.
 
     s1 = F(t, y)
-    r = norm(s1./float(max(abs(y), threshold)), Inf) + realmin() # TODO: fix type bug in max()
+    r = norm(s1./max(abs(y), threshold), Inf) + realmin() # TODO: fix type bug in max()
     h = tdir*0.8*rtol^(1/3)/r
 
     # The main loop.
@@ -79,7 +79,7 @@ function ode23(F::Function, tspan::AbstractVector, y_0::AbstractVector)
         # Estimate the error.
 
         e = h*(-5*s1 + 6*s2 + 8*s3 - 9*s4)/72
-        err = norm(e./max(float(max(abs(y), abs(ynew))), threshold), Inf) + realmin()
+        err = norm(e./max(max(abs(y), abs(ynew)), threshold), Inf) + realmin()
 
         # Accept the solution if the estimated error is less than the tolerance.
 
@@ -393,6 +393,80 @@ function ode4{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T})
     return (tspan, x)
 end
 
+#ODEROSENBROCK Solve stiff differential equations, Rosenbrock method
+#    with provided coefficients.
+function oderosenbrock{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, gamma, a, b, c)
+    h = diff(tspan)
+    x = Array(T, length(tspan), length(x0))
+    x[1,:] = x0'
+
+    solstep = 1
+    while tspan[solstep] < max(tspan)
+        ts = tspan[solstep]
+        hs = h[solstep]
+        xs = reshape(x[solstep,:], size(x0))
+        dFdx = jacobian(F, ts, xs)
+        jac = eye(size(dFdx)[1])./gamma./hs-dFdx
+
+        g = zeros(size(a)[1], length(x0))
+        g[1,:] = jac \ F(ts + b[1].*hs, xs)
+        for i = 2:size(a)[1]
+            g[i,:] = jac \ (F(ts + b[i].*hs, xs + (a[i,1:i-1]*g[1:i-1,:]).') + (c[i,1:i-1]*g[1:i-1,:]).'./hs)
+        end
+
+        x[solstep+1,:] = x[solstep,:] + b*g
+        solstep += 1
+    end
+    return(tspan, x)
+end
+
+#JACOBIAN Numerically determine Jacobian by forward finite differences
+function jacobian(F::Function, t::Number, x::AbstractVector)
+    ftx = F(t, x)
+    dFdx = zeros(length(x), length(x))
+    for j = 1:length(x)
+        dx = zeros(size(x))
+        # The 100 below is heuristic
+        dx[j] = (x[j]+(x[j]==0))./100
+        dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j]
+    end
+    return dFdx
+end
+
+# Kaps-Rentrop coefficients
+ode4s_kr(F, tspan, x0) = oderosenbrock(F,
+                                       tspan,
+                                       x0,
+                                       0.231,
+                                       [0         0        0 0
+                                        2         0        0 0
+                                        4.4524708 4.163528 0 0
+                                        4.4524708 4.163528 0 0],
+                                       [3.957037 4.624892 0.617477 1.282613],
+                                       [ 0         0         0        0
+                                        -5.071675  0         0        0
+                                         6.020153  0.159750  0        0
+                                        -1.856344 -8.505381 -2.084075 0],
+                                       )
+# Shampine coefficients
+ode4s_s(F, tspan, x0) = oderosenbrock(F,
+                                      tspan,
+                                      x0,
+                                      0.5,
+                                      [ 0    0    0 0
+                                        2    0    0 0
+                                       48/25 6/25 0 0
+                                       48/25 6/25 0 0],
+                                      [19/9 1/2 25/108 125/108],
+                                      [   0       0      0   0
+                                         -8       0      0   0
+                                        372/25   12/5    0   0
+                                       -112/125 -54/125 -2/5 0],
+                                     )
+
+# Use Shampine coefficients by default (matching Numerical Recipes)
+const ode4s = ode4s_s
+
 # ODE_MS Fixed-step, fixed-order multi-step numerical method with Adams-Bashforth-Moulton coefficients
 function ode_ms{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, order::Integer)
     h = diff(tspan)

jl/array.jl

@@ -4,9 +4,9 @@ typealias Vector{T} Array{T,1}
 typealias Matrix{T} Array{T,2}
 typealias VecOrMat{T} Union(Vector{T}, Matrix{T})
 
-typealias StridedArray{T,N}  Union(Array{T,N}, SubArray{T,N,Array{T}})
-typealias StridedVector{T} Union(Vector{T}, SubArray{T,1,Array{T}})
-typealias StridedMatrix{T} Union(Matrix{T}, SubArray{T,2,Array{T}})
+typealias StridedArray{T,N,A<:Array}  Union(Array{T,N}, SubArray{T,N,A})
+typealias StridedVector{T,A<:Array}   Union(Vector{T} , SubArray{T,1,A})
+typealias StridedMatrix{T,A<:Array}   Union(Matrix{T} , SubArray{T,2,A})
 typealias StridedVecOrMat{T} Union(StridedVector{T}, StridedMatrix{T})
 
 ## Basic functions ##

jl/inference.jl

@@ -131,7 +131,7 @@ t_func[arraysize] = (1, 2, _jl_arraysize_tfunc)
 
 function static_convert(to::ANY, from::ANY)
     if !isa(to,Tuple) || !isa(from,Tuple)
-        return (subtype(from, to) ? from : to)
+        return tintersect(from,to)
     end
     if is(to,Tuple)
         return from

jl/sparse.jl

@@ -52,15 +52,15 @@ end
 copy(S::SparseMatrixCSC) =
     SparseMatrixCSC(S.m, S.n, copy(S.colptr), copy(S.rowval), copy(S.nzval))
 
-function convert{T}(::Type{Array{T}}, S::SparseMatrixCSC{T})
+function convert{T}(::Type{Matrix{T}}, S::SparseMatrixCSC{T})
     A = zeros(T, int(S.m), int(S.n))
     for col = 1 : S.n, k = S.colptr[col] : (S.colptr[col+1]-1)
         A[S.rowval[k], col] = S.nzval[k]
     end
     return A
 end
 
-full{T}(S::SparseMatrixCSC{T}) = convert(Array{T}, S)
+full{T}(S::SparseMatrixCSC{T}) = convert(Matrix{T}, S)
 
 function sparse(A::Matrix)
     m, n = size(A)

jl/subarray.jl

@@ -304,7 +304,7 @@ strides(s::SubArray) = tuple(s.strides...)
 
 stride(s::SubArray, i::Integer) = s.strides[i]
 
-convert{T}(::Type{Ptr}, x::SubArray{T}) =
+convert{T}(::Type{Ptr{T}}, x::SubArray{T}) =
     pointer(x.parent) + (x.first_index-1)*sizeof(T)
 
 pointer(s::SubArray, i::Int) = pointer(s, ind2sub(size(s), i))

src/builtins.c

@@ -762,7 +762,7 @@ void jl_show(jl_value_t *v)
 }
 
 // comma_one prints a comma for 1 element, e.g. "(x,)"
-static void show_tuple(jl_tuple_t *t, char opn, char cls, int comma_one)
+void jl_show_tuple(jl_tuple_t *t, char opn, char cls, int comma_one)
 {
     ios_t *s = jl_current_output_stream();
     ios_putc(opn, s);
@@ -811,7 +811,7 @@ static void show_type(jl_value_t *t)
         }
         else {
             ios_write(s, "Union", 5);
-            show_tuple(((jl_uniontype_t*)t)->types, '(', ')', 0);
+            jl_show_tuple(((jl_uniontype_t*)t)->types, '(', ')', 0);
         }
     }
     else if (jl_is_seq_type(t)) {
@@ -827,7 +827,7 @@ static void show_type(jl_value_t *t)
         ios_puts(tt->name->name->name, s);
         jl_tuple_t *p = tt->parameters;
         if (p->length > 0)
-            show_tuple(p, '{', '}', 0);
+            jl_show_tuple(p, '{', '}', 0);
     }
 }
 
@@ -836,7 +836,7 @@ DLLEXPORT void jl_show_any(jl_value_t *v)
     // fallback for printing some other builtin types
     ios_t *s = jl_current_output_stream();
     if (jl_is_tuple(v)) {
-        show_tuple((jl_tuple_t*)v, '(', ')', 1);
+        jl_show_tuple((jl_tuple_t*)v, '(', ')', 1);
     }
     else if (jl_is_type(v)) {
         show_type(v);

src/gf.c

@@ -104,7 +104,7 @@ static inline int cache_match(jl_value_t **args, size_t n, jl_tuple_t *sig,
         if (jl_is_tuple(decl)) {
             // tuples don't have to match exactly, to avoid caching
             // signatures for tuples of every length
-            if (!jl_subtype(a, decl, 1))
+            if (!jl_is_tuple(a) || !jl_subtype(a, decl, 1))
                 return 0;
         }
         else if (jl_is_type_type(decl) &&
@@ -723,7 +723,17 @@ static jl_function_t *jl_mt_assoc_by_type(jl_methtable_t *mt, jl_tuple_t *tt, in
     while (m != NULL) {
         if (m->tvars!=jl_null) {
             env = jl_type_match((jl_value_t*)tt, (jl_value_t*)m->sig);
-            if (env != (jl_value_t*)jl_false) break;
+            if (env != jl_false) {
+                // parametric methods only match if all typevars are matched by
+                // non-typevars.
+                for(i=1; i < ((jl_tuple_t*)env)->length; i+=2) {
+                    if (jl_is_typevar(jl_tupleref(env,i)))
+                        break;
+                }
+                if (i >= ((jl_tuple_t*)env)->length)
+                    break;
+                env = jl_false;
+            }
         }
         else if (jl_tuple_subtype(&jl_tupleref(tt,0), nargs,
                                   &jl_tupleref(m->sig,0),
@@ -898,7 +908,8 @@ jl_methlist_t *jl_method_list_insert(jl_methlist_t **pml, jl_tuple_t *type,
     assert(jl_is_tuple(type));
     l = *pml;
     while (l != NULL) {
-        if (sigs_eq((jl_value_t*)type, (jl_value_t*)l->sig)) {
+        if (((l->tvars==jl_null) == (tvars==jl_null)) &&
+            sigs_eq((jl_value_t*)type, (jl_value_t*)l->sig)) {
             // method overwritten
             JL_SIGATOMIC_BEGIN();
             l->sig = type;
@@ -1272,6 +1283,15 @@ static void print_methlist(char *name, jl_methlist_t *ml)
     ios_t *s = jl_current_output_stream();
     while (ml != NULL) {
         ios_printf(s, "%s", name);
+        if (ml->tvars != jl_null) {
+            if (jl_is_typevar(ml->tvars)) {
+                ios_putc('{', s); jl_show((jl_value_t*)ml->tvars);
+                ios_putc('}', s);
+            }
+            else {
+                jl_show_tuple(ml->tvars, '{', '}', 0);
+            }
+        }
         jl_show((jl_value_t*)ml->sig);
         if (ml->func == NULL)  {
             // mark dummy cache entries

src/jltypes.c

@@ -1740,6 +1740,51 @@ static int jl_subtype_le(jl_value_t *a, jl_value_t *b, int ta, int morespecific,
         b==(jl_value_t*)jl_undef_type)
         return 0;
     if (!invariant && (jl_tag_type_t*)b == jl_any_type) return 1;
+
+    if (jl_is_some_tag_type(a) && jl_is_some_tag_type(b)) {
+        if ((jl_tag_type_t*)a == jl_any_type) return 0;
+        jl_tag_type_t *tta = (jl_tag_type_t*)a;
+        jl_tag_type_t *ttb = (jl_tag_type_t*)b;
+        int super=0;
+        while (tta != (jl_tag_type_t*)jl_any_type) {
+            if (tta->name == ttb->name) {
+                if (super && morespecific) {
+                    if (tta->name != jl_type_type->name)
+                        return 1;
+                }
+                if (tta->name == jl_ntuple_typename) {
+                    // NTuple must be covariant
+                    return jl_subtype_le(jl_tupleref(tta->parameters,1),
+                                         jl_tupleref(ttb->parameters,1),
+                                         0, morespecific, invariant);
+                }
+                assert(tta->parameters->length == ttb->parameters->length);
+                for(i=0; i < tta->parameters->length; i++) {
+                    jl_value_t *apara = jl_tupleref(tta->parameters,i);
+                    jl_value_t *bpara = jl_tupleref(ttb->parameters,i);
+                    if (invariant && jl_is_typevar(bpara) &&
+                        !((jl_tvar_t*)bpara)->bound) {
+                        return apara==bpara;
+                    }
+                    if (!jl_subtype_le(apara, bpara, 0, morespecific, 1))
+                        return 0;
+                }
+                return 1;
+            }
+            else if (invariant) {
+                return 0;
+            }
+            tta = tta->super; super = 1;
+        }
+        assert(!invariant);
+        if (((jl_tag_type_t*)a)->name == jl_type_type->name) {
+            // Type{T} also matches >:typeof(T)
+            if (!jl_is_typevar(jl_tparam0(a)))
+                return jl_subtype_le(jl_tparam0(a), b, 1, morespecific, 0);
+        }
+        return 0;
+    }
+
     if (jl_is_typevar(a)) {
         if (jl_is_typevar(b)) {
             return
@@ -1800,44 +1845,6 @@ static int jl_subtype_le(jl_value_t *a, jl_value_t *b, int ta, int morespecific,
     else if (jl_is_func_type(b)) {
         return 0;
     }
-
-    assert(jl_is_some_tag_type(a));
-    assert(jl_is_some_tag_type(b));
-    jl_tag_type_t *tta = (jl_tag_type_t*)a;
-    jl_tag_type_t *ttb = (jl_tag_type_t*)b;
-    int super=0;
-    while (tta != (jl_tag_type_t*)jl_any_type) {
-        if (tta->name == ttb->name) {
-            if (super && morespecific) {
-                if (tta->name != jl_type_type->name)
-                    return 1;
-            }
-            if (tta->name == jl_ntuple_typename) {
-                // NTuple must be covariant
-                return jl_subtype_le(jl_tupleref(tta->parameters,1),
-                                     jl_tupleref(ttb->parameters,1),
-                                     0, morespecific, invariant);
-            }
-            assert(tta->parameters->length == ttb->parameters->length);
-            for(i=0; i < tta->parameters->length; i++) {
-                jl_value_t *apara = jl_tupleref(tta->parameters,i);
-                jl_value_t *bpara = jl_tupleref(ttb->parameters,i);
-                if (!jl_subtype_le(apara, bpara, 0, morespecific, 1))
-                    return 0;
-            }
-            return 1;
-        }
-        else if (invariant) {
-            return 0;
-        }
-        tta = tta->super; super = 1;
-    }
-    assert(!invariant);
-    if (((jl_tag_type_t*)a)->name == jl_type_type->name) {
-        // Type{T} also matches >:typeof(T)
-        if (!jl_is_typevar(jl_tparam0(a)))
-            return jl_subtype_le(jl_tparam0(a), b, 1, morespecific, 0);
-    }
     return 0;
 }
 

src/julia.h

@@ -690,6 +690,7 @@ jl_lambda_info_t *jl_wrap_expr(jl_value_t *expr);
 
 // some useful functions
 DLLEXPORT void jl_show(jl_value_t *v);
+void jl_show_tuple(jl_tuple_t *t, char opn, char cls, int comma_one);
 
 // modules
 extern jl_module_t *jl_base_module;

test/core.jl

@@ -24,7 +24,8 @@
 @assert !isa(Array,Type{Any})
 @assert subtype(Type{ComplexPair},CompositeKind)
 @assert isa(ComplexPair,Type{ComplexPair})
-@assert subtype(Type{Ptr{None}},Type{Ptr})
+@assert !subtype(Type{Ptr{None}},Type{Ptr})
+@assert !subtype(Type{Rational{Int}}, Type{Rational})
 let T = typevar(:T)
     @assert !is(None, tintersect(Array{None},AbstractArray{T}))
     @assert  is(None, tintersect((Type{Ptr{Uint8}},Ptr{None}),
@@ -208,7 +209,7 @@ glotest()
 
 # dispatch
 begin
-    local foo, bar
+    local foo, bar, baz
     foo(x::(Any...))=0
     foo(x::(Integer...))=1
     @assert foo((:a,))==0
@@ -219,6 +220,11 @@ begin
     @assert bar((1,1,1,1)) == 1
     @assert bar((1,1,1,"a")) == 2
     @assert bar((:a,:a,:a,:a)) == 1
+
+    baz(::Type{Rational}) = 1
+    baz{T}(::Type{Rational{T}}) = 2
+    @assert baz(Rational) == 1
+    @assert baz(Rational{Int}) == 2
 end
 
 # syntax