Suggestions

src/1.JWAS/src/build_MME.jl

@@ -267,6 +267,13 @@ function getX(trm::ModelTerm,mme::MME)
     mme.mmePos  += trm.nLevels
 end
 
+function makeWeights(fw::T, useWeights::Bool, df::DataFrame, column="weights")::Vector{T} where T <: AbstractFloat
+    res = ones(typeof(fw),size(df,1))
+    if useWeights
+        res =  one(fw) ./ map((x)->(convert(typeof(fw),x)), df[:,column])  
+    end
+    return res
+end
 """
 Construct mixed model equations with
 
@@ -279,13 +286,20 @@ function getMME(mme::MME, df::DataFrame)
     if mme.mmeLhs != false
       error("Please build your model again using the function build_model().")
     end
+
+    use64 = mme.MCMCinfo == false || mme.MCMCinfo.double_precision
+
+    useWeights = mme.MCMCinfo != false && mme.MCMCinfo.heterogeneous_residuals == true
+
     #Heterogeneous residuals
-    if mme.MCMCinfo != false && mme.MCMCinfo.heterogeneous_residuals == true
-        invweights = 1 ./ convert(Array,df[!,Symbol("weights")])
+    if use64
+       invweights = makeWeights(Float64(1.0),useWeights, df, "weights")
+       mme.invweights =invweights
     else
-        invweights = ones(size(df,1))
+       invweights = makeWeights(Float32(1.0),useWeights, df, "weights")
+       mme.invweights =invweights
     end
-    mme.invweights = (mme.MCMCinfo == false || mme.MCMCinfo.double_precision ? Float64.(invweights) : Float32.(invweights))
+
     #Make incidence matrices X for each term
     for trm in mme.modelTerms
       if trm.X == false
@@ -321,7 +335,7 @@ function getMME(mme::MME, df::DataFrame)
           #rename genotype names
           mme.M[1].trait_names=mme.latent_traits
           #save omics data missing pattern
-          mme.missingPattern = .!ismissing.(convert(Matrix,df[!,mme.lhsVec]))
+          mme.missingPattern = .!ismissing.(Matrix(df[!,mme.lhsVec]))
           #replace missing data with values in yobs
           for i in mme.lhsVec      #for each omics feature
             for j in 1:size(df,1)  #for each observation
@@ -339,13 +353,14 @@ function getMME(mme::MME, df::DataFrame)
     end
     ii      = 1:length(y)
     jj      = ones(length(y))
-    vv      = ((mme.MCMCinfo == false || mme.MCMCinfo.double_precision) ? Float64.(y) : Float32.(y))
+    #vv      = ((mme.MCMCinfo == false || mme.MCMCinfo.double_precision) ? Float64.(y) : Float32.(y))
+    # This recoding should make it easy for julia to reason about the type of vv
+    vv::Union{Vector{Float64},Vector{Float32}} = use64 ? Float64.(y) : Float32.(y)
     ySparse = sparse(ii,jj,vv)
 
     #Make lhs and rhs for MME
     mme.X       = X
     mme.ySparse = ySparse
-
     if mme.nModels==1     #single-trait (lambda version)
         mme.mmeLhs = X'*Diagonal(mme.invweights)*X
         mme.mmeRhs = X'*Diagonal(mme.invweights)*ySparse
@@ -363,7 +378,6 @@ function getMME(mme::MME, df::DataFrame)
       mme.mmeLhs = X'Ri*X
       mme.mmeRhs = X'Ri*ySparse
     end
-
     #Random effects parts in MME
     if mme.nModels == 1
       #random_term.GiNew*mme.R - random_term.GiOld*mme.ROld
@@ -384,8 +398,8 @@ function getMME(mme::MME, df::DataFrame)
     #No phenotypic data for some levels of a factor in multi-trait analysis
     #e.g., y3:x3:f in https://github.com/reworkhow/JWAS.jl/blob/
     #a6f4595796b70811c0b745af525e7e0a822bb954/src/5.Datasets/data/example/phenotypes.txt
-    for i in size(mme.mmeLhs,1)
-      if mme.mmeLhs[i,i] == 0.0
+    for (i,d) in enumerate(Array(diag(mme.mmeLhs)))
+      if iszero(d)
         error("No phenotypic data for ",getNames(mme)[i])
       end
     end

src/1.JWAS/src/markers/readgenotypes.jl

@@ -103,13 +103,13 @@ function get_genotypes(file::Union{AbstractString,Array{Float64,2},Array{Float32
         #set type for each column
         ncol= length(row1)
         etv = Array{DataType}(undef,ncol)
-        fill!(etv,Float64)
+        fill!(etv,Float32)
         etv[1]=String
         close(myfile)
         #read a large genotype file
         data      = CSV.read(file,DataFrame,types=etv,delim = separator,header=false,skipto=(header==true ? 2 : 1))
-        obsID     = map(string,data[!,1])
-        genotypes = map(Float32,Matrix(data[!,2:end]))
+        obsID     = data[!,1]
+        genotypes = Matrix(data[!,2:end])
     elseif typeof(file) == DataFrames.DataFrame #Datafarme
         println("The first column in the dataframe should be individual IDs.")
         println("The data type of markers should be Number.")
@@ -119,7 +119,7 @@ function get_genotypes(file::Union{AbstractString,Array{Float64,2},Array{Float32
             markerID = string.(1:(size(file,2)-1))
         end
         obsID     = map(string,file[!,1])
-        genotypes = map(Float32,convert(Matrix,file[!,2:end]))
+        genotypes = map(Float32,Matrix(file[!,2:end]))
     elseif typeof(file) <: Union{Array{Float64,2},Array{Float32,2},Array{Any,2}} #Array
         if length(header) != (size(file,2)+1)
             header = ["id"; string.(1:size(file,2))]

src/1.JWAS/src/residual.jl

@@ -22,7 +22,7 @@ function mkRi(mme::MME,df::DataFrame,Rinv)
     nObs   = size(tstMsng,1)
     ii = Array{Int64}(undef,nObs*ntrait^2)
     jj = Array{Int64}(undef,nObs*ntrait^2)
-    vv = Array{AbstractFloat}(undef,nObs*ntrait^2)
+    vv = Array{eltype(Rinv)}(undef,nObs*ntrait^2)
     pos = 1
     for i=1:nObs
         sel = tstMsng[i,:]