ENH: special object for decision rules. #126

albop · 2017-12-20T19:13:40Z

@sglyon : here is an attempt to resolve issue #125 .
One can do:

drs = DRStore([rand(32) for i=1:3])
drs.data :: Tuple{Vector{Float64}}
drs.flat :: Vector{Float64}
maximum(abs, cat(1, drs.data...)- drs.flat) # -> 0
drs.flat[10] = 1.2
maximum(abs, cat(1, drs.data...)- drs.flat) # -> 0

The fields drs.data can contain the decision rule as we store it right now (almost, it is a tuple instead of a vector so that it is immutable) and drs.flat represents it in all states as one single vector. Trick is both fields refer to the same part of memory and are always synchronized. I believe there is no memory cleaning problem but that would require careful thinking.

I have hesitated as to how drs itself should behave. In the end, I implemented the behaviour from before (so drs[1] returns a vector) but didn't use it. The reason is to backward compatibility, but I ended up not using drs directly, using instead drs.data most of the time.
The alternative would drs to behave as drs.flat, at which case it could derive from AbstractVector. I can see how we would get some operations for free, but haven't found one very good reason. An idea would be to do drs[i] for flat indexing and drs[G(i)] or smthg like that for group indexing.

So far I have only applied it to time_iteration_direct (memory 350->20) and would like some early feedback before generalizing it. Any thought ?

sglyon

This is exciting!

I have a few concerns and questions that I think we should address before merging, but it is clear that moving in a direction like this can definitely lead to large memory efficiency gains.

sglyon · 2017-12-21T02:04:52Z

src/drstore.jl

+function DRStore(::Type{T}, dims::AbstractVector{Int}) where T
+    L = sum(dims)
+    fD = zeros(T, L)
+    ptr = pointer(fD)


Yikes 😨 ! I don't think I'v ever seen pointer and unsafe_wrap used out in in the pure Julia wild (e.g. not interfacing with C libs)

Are we confident that this is safe and works correctly? I would guess that it does given your experiments, but it makes me a little scared.

What type of penalty would we get if we were to construct D using SubArrays (e.g. using view) instead of unsafe_wrap and pointer offsets?

I think problem is not whether it is safe:

we control object construction, and the object is immutable so user cannot alter its behaviour. Of course this assumes we add consistency checks.

all arrays are linked to the same object, so one "view" cannot be destroyed without the other being destroyed. By construction, the memory we access should be freed for all arrays/subarrays at the same time.
One potential problem is actually performance. From my reptilian understanding of the discussion here (Get rid of reinterpret in its current form JuliaLang/julia#22849), the compiler can optimize further if we guarantee non-aliasing. Se we might slow down things a little bit with this approach (but not more than our current use of reinterpret).

In principle, SubArray should work, and we should be able to pass them to all routines (including interpolation routines). Right now it doesn't:

our interpolation routine assume that input are contiguous arrays (i.e. take Array as arguments). We rely on this to reinterpret arrays in both directions.

we cannot just change all arguments to AbstractArray. That would be inapropriate, because it doesn't give us any guarantee that we can reinterpret them.

SubArray have a flag which tells whether they are contiguous, but reinterpret doesn't work on them.
Intuitively, the condition for us to be able to reinterpret an array, should not for it to be contiguous but for all consecutive chunks to be a multiple of the size of all datatype we wish to reinterpret it into. I think this is the direction, which is followed in Julia 0.7, where reinterpret will return a special type of Array. It would be natural that some SubArrays are reinterpretable too, but I don't know whether this is the case yet.

Right now, I don't totally understand how the new reintepret will work in 0.7 and it certainly won't be backward compatible.
On the other hand, the current behaviour is arguably stupid, but certainly forward compatible, unless they remove pointer, which I'm sure they won't. Maybe there is a way to actually tell the compiler we are not aliasing anything. Overall there are less factors for us to take into account with that simple solution so I'm quite inclined to stay with it for now.

sglyon · 2017-12-21T02:05:55Z

src/drstore.jl

+    ptr = pointer(fD)
+    sz = sizeof(T)
+    offsets = [[0]; cumsum(dims)[1:end-1]]*sz
+    D = [unsafe_wrap(Array{T}, ptr+offsets[i], dims[i]) for i=1:length(dims)]


I think we should probably be testing that fD is the correct length, right?

unsafe_wrap will let us have elements of D defined over memory that extends beyond fD:

julia> fD = rand(4) 4-element Array{Float64,1}: 0.36243 0.0790668 0.232033 0.975159 julia> X = unsafe_wrap(Array{Float64}, pointer(fD), 5) 5-element Array{Float64,1}: 0.36243 0.0790668 0.232033 0.975159 9.88131e-324 julia> fD[:] = 1 1 julia> X 5-element Array{Float64,1}: 1.0 1.0 1.0 1.0 9.88131e-324

(they did warn us it was unsafe)

I don't see anything wrong with the value 9.88131e-324. I quite like it actually.
(joke aside: sure we need to check dimensions).

sglyon · 2017-12-21T02:07:47Z

src/drstore.jl

+
+DRStore! = DRStore
+
+function Main.copy(ds::DRStore{T}) where T


This should be Base.copy

Sure. I followed the error message without thinking...

sglyon · 2017-12-21T02:08:09Z

src/drstore.jl

+    return ds1
+end
+
+Main.copy!(ds0::DRStore{T}, ds1::DRStore{T}) where T = copy!(ds0.flat, ds1.flat)


Base.copy!

sglyon · 2017-12-21T02:08:15Z

src/drstore.jl

+Main.copy!(ds0::DRStore{T}, ds1::DRStore{T}) where T = copy!(ds0.flat, ds1.flat)
+
+
+function Main.clamp!(ds::DRStore{T}, ds_lb::DRStore{T}, ds_ub::DRStore{T}) where T


Base.clamp!

sglyon · 2017-12-21T02:11:42Z

src/algos/time_iteration_direct.jl


    while it<maxit && err>tol_η

        it += 1

        tic()

-        E_f = [zeros(Point{n_h},N) for i=1:number_of_smooth_drs(dprocess)]
+        set_values!(dr, [ds0.data...])


Is there any kind of penalty from the [ds0.data...] that appears throughout? (e.g. performance, memory, or type inference penalty)

This might be pushing things too far (and please tell me if it is), but would it make sense to have drs.data be NTuple{N,Vector{T}} so we know how many elements are in it and can use @generated functions to unroll these [drs.data...] calls?

For now I don't think we have any penalty. The elements of [ds0.data...] should be correctly inferred. set_values! itself could be improved to take a tuple as argument. I doubt unrolling would provide much benefit, as the amount of code within the loop is quite big.

albop · 2017-12-21T12:28:59Z

Thanks for the comments. So, this seems like the right direction. I'll make two changes:

rename drs.data to drs.grouped so that the object is more neutral as to which is its natural layout.
remove the indexing functions
Then I'll wait to see whether I'm hundred convinced to derive from AbstractArray{T}.
Any idea for a good name of a function to access drs.grouped[i] ? συβ(drs,i) ?

sglyon · 2017-12-21T14:21:58Z

What if we did drs[i, :] for drs.grouped[i] (or the colon in the first position -- whichever feels more natural to you)

albop · 2017-12-21T15:52:48Z

That would be a bit like a non square 2d array. Sounds very good to me !

…

On Thu, 21 Dec 2017, 14:21 Spencer Lyon, ***@***.***> wrote: What if we did drs[i, :] for drs.grouped[i] (or the colon in the first position -- whichever feels more natural to you) — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#126 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAQ5KVvx5JYIJHG33uZjBqqlZm3rbnQjks5tCmmGgaJpZM4RI1m_> .

albop · 2017-12-21T16:05:00Z

As for orientation, dolo notations are (i,n) where i is exogenous index and n endogenous one so your orientation makes sense to me. It's not classical Julia ordering but I can live with that.

…

On Thu, 21 Dec 2017, 15:52 Pablo Winant, ***@***.***> wrote: That would be a bit like a non square 2d array. Sounds very good to me ! On Thu, 21 Dec 2017, 14:21 Spencer Lyon, ***@***.***> wrote: > What if we did drs[i, :] for drs.grouped[i] (or the colon in the first > position -- whichever feels more natural to you) > > — > You are receiving this because you authored the thread. > Reply to this email directly, view it on GitHub > <#126 (comment)>, > or mute the thread > <https://github.com/notifications/unsubscribe-auth/AAQ5KVvx5JYIJHG33uZjBqqlZm3rbnQjks5tCmmGgaJpZM4RI1m_> > . >

sglyon · 2017-12-21T19:36:53Z

Ok sounds great!

albop · 2017-12-21T20:36:02Z

Quick question : what is a nice simple way to do a[:] += b[:] without allocating ?

sglyon · 2017-12-21T21:13:19Z

I believe that a .+= b would work

You can always use the trusty AXPY Blas routine:

julia> a = [1, 2, 3.0]
b =3-element Array{Float64,1}:
  1.0
 2.0
 3.0

julia> b = [10, 10, 10.0]
3-element Array{Float64,1}:
 10.0
 10.0
 10.0

julia> Base.axpy!(1, b, a)
3-element Array{Float64,1}:
 11.0
 12.0
 13.0

julia> a
3-element Array{Float64,1}:
 11.0
 12.0
 13.0

It looks like the broadcasting version is faster:

julia> function f1(a, b)
       a .+= b
       end
f1 (generic function with 1 method)

julia> function f2(a, b)
       Base.axpy!(1, b, a)
       end
f2 (generic function with 1 method)

julia> a = rand(1000);

julia> b = rand(1000);

julia> f1(a, b); @allocated f1(a, b)
96

julia> f2(a, b); @allocated f2(a, b)
0

julia> using BenchmarkTools

julia> @benchmark f1(a, b)
BenchmarkTools.Trial:
  memory estimate:  96 bytes
  allocs estimate:  4
  --------------
  minimum time:     341.165 ns (0.00% GC)
  median time:      353.771 ns (0.00% GC)
  mean time:        394.723 ns (2.26% GC)
  maximum time:     13.852 μs (94.29% GC)
  --------------
  samples:          10000
  evals/sample:     218

julia> @benchmark f2(a, b)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.381 μs (0.00% GC)
  median time:      1.403 μs (0.00% GC)
  mean time:        1.541 μs (0.00% GC)
  maximum time:     74.430 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     10

albop · 2017-12-21T21:29:24Z

Are you sure about the former ? Any ref ?

…

On Thu, 21 Dec 2017, 22:13 Spencer Lyon, ***@***.***> wrote: I believe that a .+= b would work You can always use the trusty AXPY Blas routine: julia> a = [1, 2, 3.0] b =3-element Array{Float64,1}: 1.0 2.0 3.0 julia> b = [10, 10, 10.0]3-element Array{Float64,1}: 10.0 10.0 10.0 julia> Base.axpy!(1, b, a)3-element Array{Float64,1}: 11.0 12.0 13.0 julia> a3-element Array{Float64,1}: 11.0 12.0 13.0 — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#126 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAQ5Ka-VqvbMrCa5Mk3sWlH3aPn2VkQ7ks5tCsnvgaJpZM4RI1m_> .

sglyon · 2017-12-21T22:00:56Z

Yeah I’m sure about it.

a .+= b is syntactic sugar for broadcast!(+, a, a, b) , which says (in ordrer arguments appear in the function call) broadcast the + function store the result in a and apply + to arguments a and b

julia> expand(:(a .+= b))
:((Base.broadcast!)(+, a, a, b))

albop · 2020-05-10T18:46:25Z

Getting back to it, not understanding much. Looks like my former self was a better version of me.

albop · 2022-01-04T09:59:29Z

This is now redundant, as we are using (hopefully fast) views.

ENH: special object for decision rules.

c3fa102

sglyon suggested changes Dec 21, 2017

View reviewed changes

Merge branch 'master' into albop/inplace

c5f4a3c

albop closed this Jan 4, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

ENH: special object for decision rules. #126

ENH: special object for decision rules. #126

albop commented Dec 20, 2017

sglyon left a comment

sglyon Dec 21, 2017

albop Dec 21, 2017 •

edited

Loading

sglyon Dec 21, 2017

albop Dec 21, 2017

sglyon Dec 21, 2017

albop Dec 21, 2017

sglyon Dec 21, 2017

sglyon Dec 21, 2017

sglyon Dec 21, 2017

albop Dec 21, 2017

albop commented Dec 21, 2017

sglyon commented Dec 21, 2017

albop commented Dec 21, 2017 via email

albop commented Dec 21, 2017 via email

sglyon commented Dec 21, 2017

albop commented Dec 21, 2017

sglyon commented Dec 21, 2017 •

edited

Loading

albop commented Dec 21, 2017 via email

sglyon commented Dec 21, 2017 •

edited

Loading

albop commented May 10, 2020

albop commented Jan 4, 2022


		DRStore! = DRStore

		function Main.copy(ds::DRStore{T}) where T

		Main.copy!(ds0::DRStore{T}, ds1::DRStore{T}) where T = copy!(ds0.flat, ds1.flat)


		function Main.clamp!(ds::DRStore{T}, ds_lb::DRStore{T}, ds_ub::DRStore{T}) where T

ENH: special object for decision rules. #126

ENH: special object for decision rules. #126

Conversation

albop commented Dec 20, 2017

sglyon left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

albop Dec 21, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

albop commented Dec 21, 2017

sglyon commented Dec 21, 2017

albop commented Dec 21, 2017 via email

albop commented Dec 21, 2017 via email

sglyon commented Dec 21, 2017

albop commented Dec 21, 2017

sglyon commented Dec 21, 2017 • edited Loading

albop commented Dec 21, 2017 via email

sglyon commented Dec 21, 2017 • edited Loading

albop commented May 10, 2020

albop commented Jan 4, 2022

albop Dec 21, 2017 •

edited

Loading

sglyon commented Dec 21, 2017 •

edited

Loading

sglyon commented Dec 21, 2017 •

edited

Loading