Dev (#7)

* callback * linear algebra * memory optimization * isposdef * memory optimization * minor fix, code cosmitics * svd -> svdvals * LCL' gradient optimizations * lclgf deprecate * optimizations, drop iterim data, reml function unification * deprecation, cosmetics * result structure, drop typeiii computing by default * update docs, remove structure description (not included in public API) * public API fix, deprecated functions removed * v1.0.6 * inv cholesky * update * update * update * update * update * update * update * initvar optim * simulation update * minor optimization * rbe settings init * optim * update * generalized simulation * update * update * update * update * update * cosmitics * MvNormal fix * update * update * update * fix Co-authored-by: PharmCat <[email protected]>
PharmCat · Mar 25, 2020 · b5253b8 · b5253b8 · PharmCat · Mar 25, 2020
1 parent c5922cd
commit b5253b8
Show file tree

Hide file tree

Showing 14 changed files with 476 additions and 145 deletions.
diff --git a/.travis.yml b/.travis.yml
@@ -9,14 +9,15 @@ julia:
   - 1.0
   - 1.2
   - 1.3
+  - 1.4
 
 branches:
   only:
   - master
 
 matrix:
   allow_failures:
-  - julia: 1.3
+  - julia: 1.4
 
 notifications:
   email: false

diff --git a/Project.toml b/Project.toml
@@ -3,7 +3,7 @@ name = "ReplicateBE"
 uuid = "671d9d50-c343-11e9-1a9c-fdd992682823"
 keywords = ["bioequivalence", "mixedmodel"]
 desc = "Mixed model solution for replicate designed bioequivalence study."
-version = "1.0.6"
+version = "1.0.7"
 
 [deps]
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -19,6 +19,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [extras]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -30,12 +31,13 @@ test = ["CSV", "Test", "StatsBase"]
 
 [compat]
 CategoricalArrays = "0.7"
-julia = "1.0, 1.1, 1.2"
+julia = "1.0, 1.1, 1.2, 1.3, 1.4"
 StatsBase = "0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32"
-DataFrames = "0.19"
+DataFrames = "0.19, 0.20"
 StatsModels = "0.6"
-Distributions = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21"
+Distributions = "0.20, 0.21, 0.22, 0.23"
 PDMats = "0.9"
 ForwardDiff = "0.10"
-Optim = "0.19"
+Optim = "0.19, 0.20"
 LineSearches = "7.0"
+StaticArrays = "0.11, 0.12"
diff --git a/README.md b/README.md
@@ -11,6 +11,7 @@ Mixed model solution for replicate designed bioequivalence study. This can be us
 [![codecov](https://codecov.io/gh/PharmCat/ReplicateBE.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/PharmCat/ReplicateBE.jl)
 [![Latest docs](https://img.shields.io/badge/docs-latest-blue.svg)](https://pharmcat.github.io/ReplicateBE.jl/latest/)
 [![doi](https://img.shields.io/badge/doi-10.13140%2FRG.2.2.27418.39363-blue)](https://doi.org/10.13140/RG.2.2.27418.39363)
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3543823.svg)](https://doi.org/10.5281/zenodo.3543823)
 
 Install:
 ```

diff --git a/chagelog.md b/chagelog.md
@@ -1,3 +1,15 @@
+
+- v1.0.7
+    * settings: initial 
+    * randrbeds "light" generation
+    * optimizations
+    * Distributions bump
+    * "generalized" simulation
+    * Optim 0.20
+    * Julia 1.3
+
+
+
 - v1.0.6
     * many optimizations
     * linear algebra custom functions

diff --git a/examples/simulation.jl b/examples/simulation.jl
@@ -1,6 +1,6 @@
 using ReplicateBE
 
-task = ReplicateBE.RandRBEDS(;n=24,
+task = ReplicateBE.randrbetask(;n=24,
 sequence=[1,2], design = ["T" "R" "T" "R"; "R" "T" "R" "T"],
 inter=[0.05, 0.04, 0.6], intra=[0.02, 0.04],
 intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0],
@@ -15,7 +15,7 @@ Number: 100000
 Result: 0.05078660225829358
 =#
 
-task = ReplicateBE.RandRBEDS(;n=24,
+task = ReplicateBE.randrbetask(;n=24,
 sequence=[1,2], design = ["T" "R" "T"; "R" "T" "R"],
 inter=[0.05, 0.04, 0.4], intra=[0.02, 0.04],
 intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0],
@@ -30,11 +30,53 @@ Result: 0.04997148203368122
 #E63
 =#
 
-task = ReplicateBE.RandRBEDS(;n=24,
+task = ReplicateBE.randrbetask(;n=24,
 sequence=[1,2], design = ["T" "R" "T" "R"; "R" "T" "R" "T"],
 inter=[0.05, 0.04, 0.3], intra=[0.02, 0.04],
 intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0],
 formcoef = [0.0, log(0.8)], seed = 0, dropobs = 10)
 
 result =  ReplicateBE.simulation(task; num = 2000, seed = 6564683774737, verbose = true)
 
+
+task = ReplicateBE.randrbetask(;n=24,
+sequence=[1,2], design = ["T" "R" "T" "R"; "R" "T" "R" "T"],
+inter=[0.01, 0.01, 0.9], intra=[0.09, 0.09],
+intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0],
+formcoef = [log(0.95), 0.0], seed = 0)
+
+rds    =  ReplicateBE.randrbeds(task)
+be = ReplicateBE.rbe!(rds, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence)
+
+result =  ReplicateBE.simulation(task; num = 1002, verbose = true, rsabe = false)
+
+
+
+task = ReplicateBE.randrbetask(;n=24,
+sequence=[1,2], design = ["T" "R" "T" "R"; "R" "T" "R" "T"],
+inter=2000.0, intra=[0.09, 0.09],
+intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0],
+formcoef = [log(0.95), 0.0], seed = 0)
+
+rds    =  ReplicateBE.randrbeds(task)
+result =  ReplicateBE.simulation(task; num = 1002, verbose = true, rsabe = true)
+
+
+
+
+task = ReplicateBE.randrbetask(;n=24,
+sequence=[1,1], design = ["T" "R" "T" "R"; "R" "T" "R" "T"],
+inter=[0.04, 0.04, 1.0], intra=[0.02, 0.02],
+intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0],
+formcoef = [0.0, log(0.8)], seed = 0, dropobs = 0)
+
+out = zeros(Float64, 0)
+function simfunc!(out, be)
+    push!(out, ReplicateBE.theta(be)[5])
+end
+result =  ReplicateBE.simulation!(task, out, simfunc!; num = 10, seed = 10, verbose = false)
+m      = mean(result)
+
+using Plots
+
+histogram(result)
diff --git a/src/ReplicateBE.jl b/src/ReplicateBE.jl
@@ -5,7 +5,7 @@
 
 module ReplicateBE
 
-using DataFrames, Distributions, StatsModels, StatsBase, ForwardDiff, LinearAlgebra, Random, PDMats, Optim, LineSearches, CategoricalArrays, Printf
+using DataFrames, Distributions, StatsModels, StatsBase, ForwardDiff, LinearAlgebra, Random, PDMats, Optim, LineSearches, CategoricalArrays, Printf, StaticArrays
 
     export rbe, rbe!, reml2, nobs, coef, stderror, dof, coefnum, fixed, theta, typeiii, design, show, confint, contrast, estimate, optstat, randrbeds, randrbetask
     import Base.show
@@ -16,6 +16,7 @@ using DataFrames, Distributions, StatsModels, StatsBase, ForwardDiff, LinearAlge
 const LOG2PI = log(2π)
 const MEMOPT = true
 
+include("rbesettings.jl")
 include("rbetable.jl")
 include("memalloc.jl")
 include("design.jl")

diff --git a/src/algebra.jl b/src/algebra.jl
@@ -57,7 +57,6 @@ function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::Vector, c)
             end
             @inbounds β[i] += c[n] * C[n]
         end
-
         for n = 1:p
             for m = 1:q
                 @inbounds θ[i, n] += A[m, n] * c[m]
@@ -77,7 +76,7 @@ function mulall(A::Vector, B::AbstractMatrix)
         for n = 1:q
             for m = 1:q
                 #@inbounds c[n] += B[m, n] * A[m]
-                θ += B[m, n] * A[m] * A[n]
+                @inbounds θ += B[m, n] * A[m] * A[n]
             end
             #@inbounds θ += A[n] * c[n]
         end
@@ -101,7 +100,7 @@ function mulθ₃(y::Vector, X::AbstractMatrix, β::Vector, V::AbstractMatrix, c
     end
     for n = 1:q
         for m = 1:q
-            θ += V[m, n] * (y[m] - c[m]) * (y[n] - c[n])
+            @inbounds θ += V[m, n] * (y[m] - c[m]) * (y[n] - c[n])
         end
     end
     return θ
@@ -120,12 +119,12 @@ function mulαβαt(A::AbstractMatrix, B::AbstractMatrix)
         c .= 0
         for n = 1:q
             for m = 1:q
-                c[n] +=  A[i, m] * B[n, m]
+                @inbounds c[n] +=  A[i, m] * B[n, m]
             end
         end
         for n = 1:p
             for m = 1:q
-                 mx[i, n] += A[n, m] * c[m]
+                 @inbounds mx[i, n] += A[n, m] * c[m]
             end
         end
     end
@@ -143,16 +142,17 @@ function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix)
         c .= 0
         for n = 1:q
             for m = 1:q
-                c[n] +=  A[i, m] * B[n, m]
+                @inbounds c[n] +=  A[i, m] * B[n, m]
             end
         end
         for n = 1:p
             for m = 1:q
-                 mx[i, n] += A[n, m] * c[m]
+                 @inbounds mx[i, n] += A[n, m] * c[m]
             end
         end
     end
-    mx .+ C
+    mx .+= C
+
 end
 function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, c::AbstractVector)
     q  = size(B, 1)
@@ -163,21 +163,18 @@ function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, c:
         c .= 0
         for n = 1:q
             for m = 1:q
-                try
-                c[n] +=  A[i, m] * B[n, m]
-                catch
-                    println("$q, $p, $i, $n, $m" )
-                    error()
-                end
+                @inbounds c[n] +=  A[i, m] * B[n, m]
             end
         end
         for n = 1:p
             for m = 1:q
-                 mx[i, n] += A[n, m] * c[m]
+                 @inbounds mx[i, n] += A[n, m] * c[m]
             end
         end
     end
-    mx .+ C
+    mx .+= C
+    #SMatrix{p,p,eltype(mx)}(mx)
+
 end
 """
 A * B * A' + C -> O
@@ -191,12 +188,12 @@ function mulαβαtc!(O::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C
         c .= 0
         for n = 1:q
             for m = 1:q
-                c[n] +=  A[i, m] * B[n, m]
+                @inbounds c[n] +=  A[i, m] * B[n, m]
             end
         end
         for n = 1:p
             for m = 1:q
-                 O[i, n] += A[n, m] * c[m]
+                 @inbounds O[i, n] += A[n, m] * c[m]
             end
         end
     end
@@ -211,14 +208,42 @@ function mulαβαtc!(O::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C
         c .= 0
         for n = 1:q
             for m = 1:q
-                c[n] +=  A[i, m] * B[n, m]
+                @inbounds c[n] +=  A[i, m] * B[n, m]
             end
         end
         for n = 1:p
             for m = 1:q
-                 O[i, n] += A[n, m] * c[m]
+                 @inbounds O[i, n] += A[n, m] * c[m]
             end
         end
     end
     O .+= C
 end
+
+function invchol(M)
+    q  = size(M, 1)
+    v  = zeros(eltype(M), q, q)
+    if q == 1
+        v[1,1] = 1/M[1,1]
+        return v
+    end
+    il = inv(cholesky(M).U)
+    for n = 1:q
+        for m = n:q
+            @inbounds v[n, n] += il[n, m]^2
+        end
+    end
+    for n = 1:(q-1)
+        for m = (n+1):q
+            for i = m:q
+                @inbounds v[n, m] += il[n, i] * il[m, i]
+            end
+        end
+    end
+    for n = 1:q-1
+        for m = 2:q
+            @inbounds v[m, n] = v[n, m]
+        end
+    end
+    return v
+end
diff --git a/src/deprecated.jl b/src/deprecated.jl
@@ -1,4 +1,16 @@
 
+"""
+    Feel M with zero feeled matrices
+"""
+function matvecz!(M, Zv)
+    n  = length(Zv)
+    for i = 1:n
+        M[i] = zeros(size(Zv[i])[1], size(Zv[i])[1])
+    end
+    return
+end
+
+
 struct MemAlloc{T <: AbstractFloat}
     mem1::Vector{Matrix{T}}
     mem2::Vector{Matrix{T}}