fix vcaux.res to be a matrix

OpenMendel · Jul 21, 2016 · 3a38c32 · 3a38c32
1 parent 091dc6a
commit 3a38c32
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 15 deletions.
diff --git a/src/VarianceComponentModels.jl b/src/VarianceComponentModels.jl
@@ -286,34 +286,33 @@ Auxillary data for variance component model. This can be used as pre-allocated
 intermediate variable in iterative algorithm.
 """
 immutable VarianceComponentAuxData{
- T1 <: AbstractVecOrMat,
- T2 <: AbstractMatrix,
- T3 <: AbstractVector
+ T1 <: AbstractMatrix,
+ T2 <: AbstractVector
  }
 
  res::T1 # same shape as response, p-by-d
- Xwork::T2 # nd-by-pd
- ywork::T3 # nd
- obswt::T3 # nd
+ Xwork::T1 # nd-by-pd
+ ywork::T2 # nd
+ obswt::T2 # nd
 end
 
 function VarianceComponentAuxData(vcobs::VarianceComponentVariate)
  T = eltype(vcobs)
- res = similar(vcobs.Y)
+ res = zeros(T, size(vcobs.Y, 1), size(vcobs.Y, 2))
  Xwork = zeros(T, length(vcobs.Y), nmeanparams(vcobs))
  ywork = zeros(T, length(vcobs.Y))
  obswt = zeros(T, length(vcobs.Y))
- VarianceComponentAuxData{typeof(res), typeof(Xwork), typeof(ywork)}(res,
+ VarianceComponentAuxData{typeof(Xwork), typeof(ywork)}(res,
  Xwork, ywork, obswt)
 end
 
 function VarianceComponentAuxData(vcobsrot::TwoVarCompVariateRotate)
  T = eltype(vcobsrot)
- res = similar(vcobsrot.Yrot)
+ res = zeros(T, size(vcobsrot.Yrot, 1), size(vcobsrot.Yrot, 2))
  Xwork = zeros(T, length(vcobsrot.Yrot), nmeanparams(vcobsrot))
  ywork = zeros(T, length(vcobsrot.Yrot))
  obswt = zeros(T, length(vcobsrot.Yrot))
- VarianceComponentAuxData{typeof(res), typeof(Xwork), typeof(ywork)}(res,
+ VarianceComponentAuxData{typeof(Xwork), typeof(ywork)}(res,
  Xwork, ywork, obswt)
 end
 

diff --git a/src/two_variance_component.jl b/src/two_variance_component.jl
@@ -123,7 +123,7 @@ function gradient!{T <: AbstractFloat}(
  end
  #A_mul_Bt!(N, N, vcmrot.eigvec), A_mul_B!(N, vcmrot.eigvec, N)
  N = vcmrot.eigvec * N * vcmrot.eigvec'
- copy!(sub(∇, (d^2 + 1):2d^2), N)
+ copy!(∇, d^2 + 1, N, 1, d^2)
  # gradient wrt Σ[1]
  @inbounds for j in 1:d
  λj = vcmrot.eigval[j]
@@ -140,9 +140,9 @@ function gradient!{T <: AbstractFloat}(
  end
  #A_mul_Bt!(N, N, vcmrot.eigvec), A_mul_B!(N, vcmrot.eigvec, N)
  N = vcmrot.eigvec * N * vcmrot.eigvec'
- copy!(sub(∇, 1:d^2), N)
+ copy!(∇, N)
  # scale by 0.5
- scale!(∇, convert(T, 0.5))
+ scale!(∇, 1//2)
 end # function gradient!
 
 """
@@ -1005,7 +1005,6 @@ function mm_update_Σ!{
  scale!(storage.vectors, storage.values)
  scale!(oneT ./ b1, storage.vectors)
  At_mul_B!(vcm.Σ[1], Φinv, storage.vectors)
- #A_mul_Bt!(vcm.Σ[1], vcm.Σ[1], vcm.Σ[1])
  copy!(vcm.Σ[1], A_mul_Bt(vcm.Σ[1], vcm.Σ[1]))
  # update Σ2
  scale!(b2, A2), scale!(A2, b2)
@@ -1016,7 +1015,6 @@ function mm_update_Σ!{
  scale!(storage.vectors, storage.values)
  scale!(oneT ./ b2, storage.vectors)
  At_mul_B!(vcm.Σ[2], Φinv, storage.vectors)
- #A_mul_Bt!(vcm.Σ[2], vcm.Σ[2], vcm.Σ[2])
  copy!(vcm.Σ[2], A_mul_Bt(vcm.Σ[2], vcm.Σ[2]))
 end