Access to undefined reference in unsafe_copy!

[Tony Kelman]

I'm trying to do something which is maybe a little weird, but with Julia being homoiconic it should be more feasible here than in other languages. I want a type that holds symbolic placeholders for unevaluated variables, but with mathematical operators defined to build up an expression tree around those unevaluated variables. I found an example of this in the Sims.jl package, it has an Unknown() type that does exactly what I want in the simple scalar case. If anyone can point me towards other implementations of similar functionality floating around out there, that would be cool too.

The trouble comes when I start trying to use these Unknown() objects as the nonzero elements in a sparse matrix. Here's my first crack at that:

julia> using Sims

... some warnings

julia> A = sparse(1:2, 1:2, [Unknown() for i=1:2])

ERROR: no method convert(Type{Unknown{DefaultUnknown}}, Int64) in sparse at sparse/csparse.jl:30

Okay, zero(Tv) doesn't work for the Unknown type (I have a feeling the backtrace should be longer?). How about

julia> A = sparse(1:2, 1:2, Any[Unknown() for i=1:2])

ERROR: type: non-boolean (MExpr) used in boolean context in sparse at sparse/csparse.jl:30

Same line, different problem. The != operator produces unevaluated expression objects instead of boolean, not too surprising. I can overwrite that for my purposes.

julia> !=(::Unknown, ::Number) = true

julia> A = sparse(1:2, 1:2, Any[Unknown() for i=1:2])

2x2 sparse matrix with 2 Any entries:

        [1, 1]  =  <<`549`,0.0>>

        [2, 2]  =  <<`550`,0.0>>

Awesome. But will sparse matrix multiplication work with this thing?

julia> A * A

ERROR: access to undefined reference in unsafe_copy! at array.jl:41

Now I'm confused. This was on Windows from a recent binary installer,

Julia Version 0.3.0-prerelease+1570

Commit a6f6461* (2014-02-14 21:07 UTC)

Platform Info:

  System: Windows (x86_64-w64-mingw32)

  CPU: Intel(R) Core(TM) i7-2630QM CPU @ 2.00GHz

  WORD_SIZE: 64

  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY)

  LAPACK: libopenblas

  LIBM: libopenlibm

If I go over to a from-source Linux build,

Julia Version 0.3.0-prerelease+1622

Commit eeb2b00* (2014-02-18 07:07 UTC)

Platform Info:

  System: Linux (x86_64-redhat-linux6E)

  CPU: Intel(R) Xeon(R) CPU           E5410  @ 2.33GHz

  WORD_SIZE: 64

  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY)

  LAPACK: libopenblas

  LIBM: libopenlibm

then the backtrace for that last error is slightly different, it says "in * at linalg/sparse.jl:171." So, something having to do with splice! then? Anyone have any recommendations how to figure out what's going on?

Thanks,

Tony

[Mauro]

> I'm trying to do something which is maybe a little weird, but with Julia
> being homoiconic it should be more feasible here than in other languages. I
> want a type that holds symbolic placeholders for unevaluated variables, but
> with mathematical operators defined to build up an expression tree around
> those unevaluated variables.
>
> I found an example of this in the Sims.jl package, it has an Unknown()
> type that does exactly what I want in the simple scalar case. If
> anyone can point me towards other implementations of similar
> functionality floating around out there, that would be cool too.

I played around with a type which accumulates the operations applied to it
and then can be evaluated later:
https://bitbucket.org/maurow/delayednumber.jl/src/a80825a073553c28616b2d110cfbb914de216c11/DelayedAny.jl?at=master
For example:

de = DeAny() # uses :SYM
de2 = (de + 5)*7*9
evalDe(de2, 5) # now it gets evaluated with de2=5

At a cursory glance this is quite similar to Sims.jl, but I'm sure much
much less sophisticated.

Then, I had a look into your troubles: first, using the DeAny type
instead of Unknown() doesn't work either. The error is again from that
line sparse/csparse.jl:30 which calls `zero` on the datatype. Once I
extended `zero` for DeAny, I also arrive at `ERROR: type:
non-boolean (DeAny) used in boolean context`. Once I extend !=, as you
did, it works to construct a sparse array. I chased it a bit further
but got stuck soon after. Have a look at above link if interested.

It does work when doing this with a normal array, again see above link.

Anyway, is it wise to try this? Maybe you could have the Unkown one
abstraction layer higher and only build the sparse matrix once you're
evaluating it?

[Tony Kelman]

Cool, thanks for the links. So a few folks have looked at doing this kind of thing, but maybe not specifically with sparse matrices yet :)

The application I'm most interested in is sparse nonlinear constrained optimization, but the same pattern shows up lots of places: you're computing a sparse matrix with constant nonzero structure, but many iterations in a row with changing nonzero values.

What would be perfect here is symbolic assembly of the initial nonzero structure, into a fast function straight from your vector input to the vector of nonzero values in the sparse matrix, without any interpreted temporaries in between. Miles and the other JuliaOpt guys are working on doing this with reverse-mode AD, but I've got an alternate formulation where I can express the input/output relationship in terms of spmm, with one matrix constant and the second matrix with nonzeros that are simple functions of the input.

So in my case I'm going for something a little like PETSc's separation of MatMatMultSymbolic (http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMatMultSymbolic.html) from MatMatMultNumeric. I feel like in Julia, symbolic can really mean symbolic, actually populating the array of nonzeros with the real expressions.

-Tony

[Tony Kelman]

I think I figured it out. Your code gives a much nicer stacktrace than Unknown() from Sims.jl was for some reason:

ERROR: access to undefined reference

 in unsafe_copy! at array.jl:41

 in copy! at array.jl:51

 in getindex at array.jl:296

 in splice! at array.jl:643

 in * at linalg/sparse.jl:171

She's actually completely right. Sparse matmul allocates space as it goes, in reasonable-sized chunks at a time. At the end, it trims off the extra unneeded elements with splice!. The problem is with these Any type arrays we're trying to use, the extra allocated space is actually #undef until it gets written to. Since splice! returns the removed elements, that first line of splice is indeed trying to access an undefined reference.

Minimal test case, no funky symbolic expression stuff required:

julia> sparse(1:2, 1:2, Any[1;2])^2

I think the solution is to replace uses of splice! with deleteat! in linalg/sparse.jl. Testing this now and will put together a PR if it works.

-Tony

[Mauro]

My approach to delayed execution is about 10x slower than the normal
execution, this might be a problem for you. Also, usually all
operations done to matrix entries are the same. If that is the case for
you, then you don't need to have each entry a delayed datatype. Instead
have the one delayed datatype represent the whole matrix:

de = DeAny()
de2 = (de^2)*7*9
a=sparse(1:2, 1:2, [1,2])
evalDe(de2, a)