I could live with the ugly but the problem is that I would also want -!, *!, /!, ^! ...
That’s easy, you can just overload Base.:-(::A, ::Not), Base.:*(::A, ::Not), etc. Just keep in mind that this only really works if you actually own the type on the right.
1 Like
Yes, this is to be applied to my own types.
Thanks, will try it.
1 Like
If you are trying to define in-place versions of elementary operations like +, presumably because you are working with large data structures and want to avoid heap allocations, then realize that it is usually much more efficient to combine multiple simple operations in a single “fused” operation (for cache locality and other reasons), so you might want to re-think your approach.
This is precisely the reason for Julia’s fused broadcasting semantics: @. x = x * y - 3x + 2/z fuses into a single loop acting in-place on x for arrays x,y,z. It is perfectly possible to use this syntax to define “fused” operations on your own types.
6 Likes
Well, this something I would like to offer the users of my package. Broadcasting is nice but so would be !functions but I’m not being able to implement it because I can’t see how to access the elements of A & B in above example.
I’ll look at the broadcasting too.
BTW, I’m curious about the += operator. Where is it defined?
It’s not a generic operator. a += b is lowered to a = a + b, there’s no way to overload that. If you use the broadcasting machinery like @stevengj suggested, you could intercept a .= ... and therefore also a .+= ... though, through Base.materialize! .
2 Likes
The broadcast machinery looks quite complex, and to start with the type in question does not have an axes as mentioned necessary
julia> G1 = mat2grid([0.0 1; 2 3]);
julia> axes(G1)
ERROR: MethodError: no method matching size(::GMT.GMTgrid)
Perhaps I’m missing something, but your type looks to me like a special matrix type, so is there a reason why you wouldn’t want to define size on it and implement the AbstractArray interface for it?
1 Like
No, nothing against, only ignorance. But again implementing that interface is not obvious. The Squares example in the docs is too simplistic. Would I have to declare only the type as <: AbstractArray or do I have to do it also for the array members? But I want some of them to remain Float64 only.
(and thanks for the advices)
Yes, your type should subtype AbstractArray{T,N}, where T is the element type and N the number of dimensions. It doesn’t matter what types the members of your custom array type are, that’s not important for dispatch. You probably want to implement at least Base.size and Base.getindex, which should already get you on your way. You can then try doing things like broadcasting, usually you should be able to see from the stacktraces, which additional functions you still need to implement to support a certain feature, but I think there is something in the manual about that as well.
2 Likes
OK, it was not that easy but I finally managed to implement broadcasting. I have more questions but will leave it to another post. I found, however, a strange behavior that even looks like a bug. The implementation is here
This works like expected
julia> G1 = mat2grid([0.0 pi/4; pi pi/2])
2×2 GMT.GMTgrid{Float64,2}:
0.0 0.785398
3.14159 1.5708
julia> G2 = cos.(G1)
2×2 GMT.GMTgrid{Float64,2}:
1.0 0.707107
-1.0 6.12323e-17
julia> acos.(G2)
2×2 GMT.GMTgrid{Float64,2}:
0.0 0.785398
3.14159 1.5708
julia> cos(G1)
ERROR: MethodError: no method matching exp!(::GMT.GMTgrid{Complex{Float64},2})
Closest candidates are:
exp!(::StridedArray{T, 2}) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\dense.jl:554
Stacktrace:
[1] cos(::GMT.GMTgrid{Float64,2}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\dense.jl:808
[2] top-level scope at REPL[6]:1
but this, which takes quite some time to compute at first time, does not error like the above last case and gives the following result. This is very dangerous because it results from a mistake very easy to make.
julia> acos(G2)
2×2 Array{Complex{Float64},2}:
0.651095+5.55112e-17im -0.71526+5.55112e-17im
1.01153+0.0im 1.66263-2.22045e-16im
It’s giving you the correct answer — the acos of a square matrix is perfectly well defined and that is what it is giving you.
3 Likes
It’s not the correct answer. The correct is given by acos.(). See, not even the type is correct. And it’s not only acos
julia> G1 = mat2grid([0.0 pi/4; pi pi/2])
2×2 GMT.GMTgrid{Float64,2}:
0.0 0.785398
3.14159 1.5708
julia> G2 = cos.(G1)
2×2 GMT.GMTgrid{Float64,2}:
1.0 0.707107
-1.0 6.12323e-17
julia> acos.(G2)
2×2 GMT.GMTgrid{Float64,2}:
0.0 0.785398
3.14159 1.5708
julia> acos(G2)
2×2 Array{Complex{Float64},2}:
0.651095+5.55112e-17im -0.71526+5.55112e-17im
1.01153+0.0im 1.66263-2.22045e-16im
julia> using Statistics
julia> mean.(G1)
2×2 GMT.GMTgrid{Float64,2}:
0.0 0.785398
3.14159 1.5708
julia> mean(G1)
1.3744467859455345
julia> std(G1)
1.3413227186681242
I think you have a fundamental misunderstanding here of linear algebra or are using the word ‘correct’ to instead mean ‘the thing that I meant’.
What you are seeing is no different from * versus .*.
julia> A = [1 2; 3 4]; B = [5 6; 7 8];
julia> A * B
2×2 Array{Int64,2}:
19 22
43 50
julia> A .* B
2×2 Array{Int64,2}:
5 12
21 32
3 Likes
OK, I see what you mean. It’s computing the acos matrix of the type array member, though I was expecting to get the composite type on output as well.
julia> acos(G2.z)
2×2 Array{Complex{Float64},2}:
0.651095+5.55112e-17im -0.71526+5.55112e-17im
1.01153+0.0im 1.66263-2.22045e-16im
julia> acos(G2)
2×2 Array{Complex{Float64},2}:
0.651095+5.55112e-17im -0.71526+5.55112e-17im
1.01153+0.0im 1.66263-2.22045e-16im
but then why the same does not hold with cos?
julia> cos(G2.z)
2×2 Array{Float64,2}:
0.827518 -0.365429
0.516794 1.34431
julia> cos(G2)
ERROR: MethodError: no method matching GMT.GMTgrid(::String, ::String, ::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Float64, ::String, ::String, ::String, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Complex{Float64},2}, ::String, ::String, ::String, ::String)
Closest candidates are:
GMT.GMTgrid(::String, ::String, ::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Union{Float32, Float64}, ::String, ::String, ::String, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{T,N}, ::String, ::String, ::String, ::String) where {T<:Real, N} at C:\Users\joaqu\.julia\dev\GMT\src\gmt_main.jl:2
Stacktrace:
Broadcasting is most likely still the best solution here, but note that in Julia 1.6, ꜝ will be a valid operator suffix, so you can use +ꜝ, -ꜝ, *ꜝ, etc. as infix or prefix operators, without resorting to any hacks:
3 Likes