-
Notifications
You must be signed in to change notification settings - Fork 48
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
rankUpdate for BlockDiagonal by Dense #447
Conversation
Codecov Report
@@ Coverage Diff @@
## master #447 +/- ##
==========================================
+ Coverage 92.74% 92.85% +0.11%
==========================================
Files 23 23
Lines 1722 1736 +14
==========================================
+ Hits 1597 1612 +15
+ Misses 125 124 -1
Continue to review full report at Codecov.
|
@kleinschmidt While working on this, I discovered a second issue. These two formulae should give the same model: julia> fnested = @formula(Thickness ~ 1 + Source + (1+Source|Lot/Wafer))
FormulaTerm
Response:
Thickness(unknown)
Predictors:
1
Source(unknown)
(Source,Lot,Wafer)->(1 + Source) | Lot / Wafer
julia> f = @formula(Thickness ~ 1 + Source + (1+Source|Lot) + (1+Source|Lot&Wafer))
FormulaTerm
Response:
Thickness(unknown)
Predictors:
1
Source(unknown)
(Source,Lot)->(1 + Source) | Lot
(Source,Lot,Wafer)->(1 + Source) | Lot & Wafer But julia> fit(MixedModel, fnested, MixedModels.dataset(:oxide))
ERROR: MethodError: no method matching +(::Tuple{InterceptTerm{true},CategoricalTerm{DummyCoding,String,1}}, ::Tuple{RandomEffectsTerm,RandomEffectsTerm})
Closest candidates are:
+(::Any, ::Any, ::Any, ::Any...) at operators.jl:538
+(::MutableArithmetics.Zero, ::Any) at /home/XXX/.julia/packages/MutableArithmetics/0tlz5/src/rewrite.jl:52
+(::Tuple{Vararg{Union{AbstractTerm, Tuple{Vararg{AbstractTerm,N}} where N},N}} where N, ::AbstractTerm) at /home/XXX/.julia/packages/StatsModels/cpRh6/src/terms.jl:404
...
Stacktrace:
[1] +(::InterceptTerm{true}, ::CategoricalTerm{DummyCoding,String,1}, ::Tuple{RandomEffectsTerm,RandomEffectsTerm}) at ./operators.jl:538
[2] sum(::Tuple{InterceptTerm{true},CategoricalTerm{DummyCoding,String,1},Tuple{RandomEffectsTerm,RandomEffectsTerm}}) at ./tuple.jl:396
[3] apply_schema(::Tuple{ConstantTerm{Int64},Term,FunctionTerm{typeof(|),var"#107#109",(:Source, :Lot, :Wafer)}}, ::MixedModels.MultiSchema{StatsModels.FullRank}, ::Type{T} where T) at /home/XXX/Work/MixedModels.jl/src/schema.jl:30
[4] apply_schema(::FormulaTerm{Term,Tuple{ConstantTerm{Int64},Term,FunctionTerm{typeof(|),var"#107#109",(:Source, :Lot, :Wafer)}}}, ::StatsModels.Schema, ::Type{LinearMixedModel}) at /home/XXX/Work/MixedModels.jl/src/schema.jl:53
[5] LinearMixedModel(::FormulaTerm{Term,Tuple{ConstantTerm{Int64},Term,FunctionTerm{typeof(|),var"#107#109",(:Source, :Lot, :Wafer)}}}, ::NamedTuple{(:Source, :Lot, :Wafer, :Site, :Thickness),Tuple{Arrow.DictEncoded{String,UInt32,SentinelArrays.ChainedVector{String,Arrow.List{String,Int32,Array{UInt8,1}}}},Arrow.DictEncoded{String,UInt32,SentinelArrays.ChainedVector{String,Arrow.List{String,Int32,Array{UInt8,1}}}},Arrow.DictEncoded{String,UInt32,SentinelArrays.ChainedVector{String,Arrow.List{String,Int32,Array{UInt8,1}}}},Arrow.DictEncoded{String,UInt32,SentinelArrays.ChainedVector{String,Arrow.List{String,Int32,Array{UInt8,1}}}},Arrow.Primitive{Float64,Array{Float64,1}}}}; contrasts::Dict{Symbol,Any}, wts::Array{Any,1}) at /home/XXX/Work/MixedModels.jl/src/linearmixedmodel.jl:67 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
A few stylistic comments, none of which are terribly important. Do with them what you will.
src/linalg/rankUpdate.jl
Outdated
for i in 1:blksize, j in 1:blksize | ||
iind = ioffset + i | ||
jind = joffset + j | ||
Arow = view(A, iind, :) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would probably write this as an explicit loop over axes(A, 2)
and avoid the construction of views.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
So something like this:
AtAij = 0
for idx in axes(A, 2)
# because the second is actually A', we swap index orders
AtAij += A[iind,idx] * A[jind,idx]
end
Cdat[i,j,k] += α * AtAij
?
I also thought about writing this as
AtAij = sum(A[iind,idx] * A[jind,idx] for idx in axes(A,2))
but then I realized the construction of the generator is probably no better (and maybe worse) than the construction of views.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Actually what I had in mind is something like
for idx in axes(A, 2)
Cdat[i,j,k] += α * A[iind, idx] * A[jind, idx]
end
I'm pretty sure the compiler will detect that the location Cdat[i,j,k]
is a loop invariant.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah, interesting.
That might be a decent research project if we ever get an intern -- taking a look at @code_llvm
and @code_native
for some of our loops and seeing where the Julia compiler and LLVM are doing nice vectorization and other loop transformations for.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
An alternative is to resize Cdat
to be a matrix with blksize
rows. Then the columns match up with the columns of C
. Also there is an intriguing function called fldmod1
which makes it easy to step over the blocks but I have forgotten the details. I'll take a look.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If we're going to resize in that way -- would it make sense to change the internal storage used by UniformBlockDiagonal
to always store in that format?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There are times that you want the diagonal blocks to be accessed as k by k matrices. The most important part is in evaluating Lambda'Z'ZLambda + I
. Also, the solving of multiple triangular systems to produce the [i,1]
blocks of L for i > 1
. I think it is easier to resize from 3 dimensions to 2 when we need it, rather than the other way around.
Co-authored-by: Douglas Bates <[email protected]>
Co-authored-by: Douglas Bates <[email protected]>
Co-authored-by: Douglas Bates <[email protected]>
Co-authored-by: Douglas Bates <[email protected]>
That's really strange, I don't see why that wouldn't work. It looks like it's not failing at the the random slope point, but rather at the point where it's trying to combine the two RE terms with the FE terms...but that path should already be hit by the existing tests, so I'm confused. I'll try to take a look when I get a chance. |
@kleinschmidt Thanks! Glad to know at least the problem isn't something obvious I should have been able to fix in 10 minutes... |
Another one from #445.