From 53881df3b62f35f4ac873d3686c33cafb86a9ac0 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Tue, 23 Jan 2024 21:28:54 +0300 Subject: [PATCH 1/7] experimental --- src/varstruct.jl | 46 ++++++++++++++++++++++++++++++++++++---------- 1 file changed, 36 insertions(+), 10 deletions(-) diff --git a/src/varstruct.jl b/src/varstruct.jl index d87dd76..4b319cf 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -1,5 +1,21 @@ const CType = Union{FunctionTerm{typeof(+), Vector{Term}}, FunctionTerm{typeof(*), Vector{Term}}, FunctionTerm{typeof(&), Vector{Term}}} +import StatsModels: ContrastsMatrix, AbstractContrasts, modelcols + +mutable struct RawCoding <: AbstractContrasts +end +function StatsModels.ContrastsMatrix(contrasts::RawCoding, levels::AbstractVector{T}) where T + ContrastsMatrix(ones(1,1), + ["levels"], + levels, + contrasts) +end +function StatsModels.modelcols(t::CategoricalTerm{RawCoding, T, N}, d::NamedTuple) where T where N + #v = d[t.sym] + #reshape(v, length(v), 1) + d[t.sym] +end + ################################################################################ # @covstr macro ################################################################################ @@ -154,7 +170,7 @@ end """ Covarince structure. """ -struct CovStructure{T} <: AbstractCovarianceStructure +struct CovStructure{T, T2} <: AbstractCovarianceStructure # Random effects random::Vector{VarEffect} # Repearted effects @@ -180,7 +196,7 @@ struct CovStructure{T} <: AbstractCovarianceStructure # unit range z column range for each random effect zrndur::Vector{UnitRange{Int}} # repeated effect parametrization matrix - rz::Vector{Matrix{T}} + rz::Vector{Matrix{T2}} # size 2 of z/rz matrix q::Vector{Int} # total number of parameters in each effect @@ -209,7 +225,7 @@ struct CovStructure{T} <: AbstractCovarianceStructure z = Matrix{Float64}(undef, rown, 0) #subjz = Vector{BitMatrix}(undef, alleffl) dicts = Vector{Dict}(undef, alleffl) - # + # unit range z column range for each random effect zrndur = Vector{UnitRange{Int}}(undef, length(random)) # Number of random effects rn = length(random) @@ -218,7 +234,7 @@ struct CovStructure{T} <: AbstractCovarianceStructure # Number of repeated effects rpn = length(repeated) # Z Matrix for repeated effect - rz = Vector{Matrix{Float64}}(undef, rpn) + # rz = Vector{Matrix{Float64}}(undef, rpn) # #Theta parameter type ct = Vector{Symbol}(undef, 0) @@ -239,7 +255,12 @@ struct CovStructure{T} <: AbstractCovarianceStructure if length(random[i].coding) == 0 fill_coding_dict!(random[i].model, random[i].coding, data) end - schema[i] = apply_schema(random[i].model, StatsModels.schema(data, random[i].coding)) + #data_ = data[] + if isa(random[i].covtype.s, ZERO) + schema[i] = InterceptTerm{false}() + else + schema[i] = apply_schema(random[i].model, StatsModels.schema(data, random[i].coding)) + end ztemp = modelcols(MatrixTerm(schema[i]), data) q[i] = size(ztemp, 2) csp = covstrparam(random[i].covtype.s, q[i]) @@ -261,16 +282,18 @@ struct CovStructure{T} <: AbstractCovarianceStructure append!(emap, fill(i, t[i])) rtn += t[i] end - + + rz_ = Vector{Matrix}(undef, rpn) # REPEATED EFFECTS for i = 1:length(repeated) if length(repeated[i].coding) == 0 fill_coding_dict!(repeated[i].model, repeated[i].coding, data) end - + schema[rn+i] = apply_schema(repeated[i].model, StatsModels.schema(data, repeated[i].coding)) - rz[i] = modelcols(MatrixTerm(schema[rn+i]), data) + #rz_[i] = reduce(hcat, modelcols(schema[rn+i], data)) + rz_[i] = modelcols(MatrixTerm(schema[rn+i]), data) symbs = StatsModels.termvars(repeated[i].subj) if length(symbs) > 0 cdata = tabcols(data, symbs) # Tuple(Tables.getcolumn(Tables.columns(data), x) for x in symbs) @@ -281,7 +304,7 @@ struct CovStructure{T} <: AbstractCovarianceStructure end sn[rn+i] = length(dicts[rn+i]) - q[rn+i] = size(rz[i], 2) + q[rn+i] = size(rz_[i], 2) csp = covstrparam(repeated[i].covtype.s, q[rn+i]) t[rn+i] = sum(csp) tr[rn+i] = UnitRange(sum(t[1:rn+i-1]) + 1, sum(t[1:rn+i-1]) + t[rn+i]) @@ -290,6 +313,9 @@ struct CovStructure{T} <: AbstractCovarianceStructure # emap append!(emap, fill(rn+i, t[rn+i])) end + T2 = typejoin(eltype.(rz_)...) + rz = Vector{Matrix{T2}}(undef, rpn) + rz .= rz_ # Theta length tl = sum(t) ######################################################################## @@ -351,7 +377,7 @@ struct CovStructure{T} <: AbstractCovarianceStructure end esb = EffectSubjectBlock(sblock, nblock) ####################################################################### - new{eltype(z)}(random, repeated, schema, rcnames, blocks, rn, rtn, rpn, z, esb, zrndur, rz, q, t, tr, tl, ct, emap, sn, maxn) + new{eltype(z), T2}(random, repeated, schema, rcnames, blocks, rn, rtn, rpn, z, esb, zrndur, rz, q, t, tr, tl, ct, emap, sn, maxn) end end ############################################################################### From 1e294ece2f3b54e46533fb9146359f8e8161b28d Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 26 Jan 2024 02:26:19 +0300 Subject: [PATCH 2/7] row col --- Project.toml | 2 +- docs/src/api.md | 14 +++++++++++++ docs/src/custom.md | 28 ++++++++++++++++++++++++++ src/Metida.jl | 2 +- src/rmat.jl | 5 +++++ src/utils.jl | 2 +- src/varstruct.jl | 49 +++++++++++++++++++++++++++++++--------------- test/devtest.jl | 19 +++++++++++++++++- test/test.jl | 14 ++++++++++++- 9 files changed, 114 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index c91d01e..86d8102 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,7 @@ uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2" keywords = ["lenearmodel", "mixedmodel"] desc = "Mixed-effects models with flexible covariance structure." authors = ["Vladimir Arnautov "] -version = "0.15.0" +version = "0.15.1" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" diff --git a/docs/src/api.md b/docs/src/api.md index 04a525e..b2c824b 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -400,6 +400,13 @@ Metida.typeiii Metida.MILMM ``` +### Metida.RawCoding + +```@docs +Metida.RawCoding +``` + + ## Not API functions ### Metida.contrast @@ -443,3 +450,10 @@ Metida.tname ```@docs Metida.raneflenv ``` + +### Metida.edistance + +```@docs +Metida.edistance +``` + diff --git a/docs/src/custom.md b/docs/src/custom.md index 2dcfac1..72fb0a0 100644 --- a/docs/src/custom.md +++ b/docs/src/custom.md @@ -90,6 +90,7 @@ Example: ```@example lmmexample using Metida, DataFrames, CSV, CategoricalArrays +spatdf = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "spatialdata.csv"); types = [Int, Int, String, Float64, Float64, Float64, Float64, Float64]) |> DataFrame ftdf = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "1fptime.csv"); types = [String, String, Float64, Float64]) |> DataFrame df0 = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "df0.csv"); types = [String, String, String, String,Float64, Float64, Float64]) |> DataFrame @@ -143,3 +144,30 @@ repeated = Metida.VarEffect(Metida.@covstr(period|subject), CustomCovarianceStru ) Metida.fit!(lmm) ``` + +### Custom distance estimation for spatial structures + +If you want to use coordinates or some other structures for distance estimation you can define method [`Metida.edistance`](@ref) to calculate distance: + +```@example lmmexample +function Metida.edistance(mx::AbstractMatrix{<:CartesianIndex}, i::Int, j::Int) + return sqrt((mx[i, 1][1] - mx[j, 1][1])^2 + (mx[i, 1][2] - mx[j, 1][2])^2) +end +``` + +For example this method returns distance between two vectors represented as `CartesianIndex`. + +Make vector of `CartesianIndex`: + +```@example lmmexample +spatdf.ci = map(x -> CartesianIndex(x[:x], x[:y]), eachrow(spatdf)) +``` + +Then use new collumn as "raw" variable with [`Metida.RawCoding`](@ref) contrast and fit the model: + +```@example lmmexample +lmm = Metida.LMM(@formula(r2 ~ f), spatdf; + repeated = Metida.VarEffect(Metida.@covstr(ci|1), Metida.SPEXP; coding = Dict(:ci => Metida.RawCoding())), + ) +Metida.fit!(lmm) +``` \ No newline at end of file diff --git a/src/Metida.jl b/src/Metida.jl index ca1059a..70685a3 100644 --- a/src/Metida.jl +++ b/src/Metida.jl @@ -32,7 +32,7 @@ TOEPH, HeterogeneousToeplitz, TOEPHP, HeterogeneousToeplitzParameterized, SPEXP, SpatialExponential, SPPOW, SpatialPower, -SPGAU, SpatialGaussian, +SPGAU, SpatialGaussian, RawCoding, UN, Unstructured, CovarianceType, fit, fit!, LMM, VarEffect, theta, logreml, m2logreml, thetalength, dof_satter, dof_contain, rankx, caic, lcontrast, typeiii, estimate, contrast, diff --git a/src/rmat.jl b/src/rmat.jl index 15b0154..04fcaf7 100644 --- a/src/rmat.jl +++ b/src/rmat.jl @@ -167,6 +167,11 @@ end return sqrt(sum) end =# +""" + edistance(mx::AbstractMatrix{T}, i::Int, j::Int) where T + +Distance between vector mx[i, :] and mx[j, :]. +""" function edistance(mx::AbstractMatrix{T}, i::Int, j::Int) where T sum = zero(T) @inbounds for c = 1:size(mx, 2) diff --git a/src/utils.jl b/src/utils.jl index 64ca46e..2bab8fa 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -31,7 +31,7 @@ function nterms(rhs::Union{Tuple{Vararg{AbstractTerm}}, Nothing, AbstractTerm}) p end """ - Rerm name. + Term name. """ tname(t::AbstractTerm) = "$(t.sym)" tname(t::InteractionTerm) = join(tname.(t.terms), " & ") diff --git a/src/varstruct.jl b/src/varstruct.jl index 4b319cf..cf1c3ff 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -2,6 +2,11 @@ const CType = Union{FunctionTerm{typeof(+), Vector{Term}}, FunctionTerm{typeof(* import StatsModels: ContrastsMatrix, AbstractContrasts, modelcols +""" + mutable struct RawCoding <: AbstractContrasts + +Contrast for CategoricalTerm to get collumn "as it is" for model matrix. +""" mutable struct RawCoding <: AbstractContrasts end function StatsModels.ContrastsMatrix(contrasts::RawCoding, levels::AbstractVector{T}) where T @@ -255,22 +260,30 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure if length(random[i].coding) == 0 fill_coding_dict!(random[i].model, random[i].coding, data) end - #data_ = data[] + if isa(random[i].model, ConstantTerm) # if only ConstantTerm in the model - data_ - first is collumn (responce) + data_ = data[[first(keys(data))]] + else + data_ = data[StatsModels.termvars(random[i].model)] # only collumns for model + end if isa(random[i].covtype.s, ZERO) schema[i] = InterceptTerm{false}() + zsize = 0 else - schema[i] = apply_schema(random[i].model, StatsModels.schema(data, random[i].coding)) + schema[i] = apply_schema(random[i].model, StatsModels.schema(data_, random[i].coding)) + ztemp = modelcols(MatrixTerm(schema[i]), data_) + z = hcat(z, ztemp) + zsize = size(ztemp, 2) end - ztemp = modelcols(MatrixTerm(schema[i]), data) - q[i] = size(ztemp, 2) + + q[i] = zsize csp = covstrparam(random[i].covtype.s, q[i]) t[i] = sum(csp) - z = hcat(z, ztemp) + fillur!(zrndur, i, q) fillur!(tr, i, t) symbs = StatsModels.termvars(random[i].subj) if length(symbs) > 0 - cdata = tabcols(data, symbs) # Tuple(Tables.getcolumn(Tables.columns(data), x) for x in symbs) + cdata = tabcols(data, symbs) # Tuple(Tables.getcolumn(Tables.columns(data_), x) for x in symbs) dicts[i] = Dict{Tuple{eltype.(cdata)...}, Vector{Int}}() indsdict!(dicts[i], cdata) else @@ -290,26 +303,30 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure if length(repeated[i].coding) == 0 fill_coding_dict!(repeated[i].model, repeated[i].coding, data) end + if isa(repeated[i].model, ConstantTerm) # if only ConstantTerm in the model - data_ - first is collumn (responce) + data_ = data[[first(keys(data))]] + else + data_ = data[StatsModels.termvars(repeated[i].model)] # only collumns for model + end - schema[rn+i] = apply_schema(repeated[i].model, StatsModels.schema(data, repeated[i].coding)) + schema[rn + i] = apply_schema(repeated[i].model, StatsModels.schema(data_, repeated[i].coding)) #rz_[i] = reduce(hcat, modelcols(schema[rn+i], data)) - rz_[i] = modelcols(MatrixTerm(schema[rn+i]), data) + rz_[i] = modelcols(MatrixTerm(schema[rn+i]), data_) symbs = StatsModels.termvars(repeated[i].subj) if length(symbs) > 0 - cdata = tabcols(data, symbs) # Tuple(Tables.getcolumn(Tables.columns(data), x) for x in symbs) - dicts[rn+i] = Dict{Tuple{eltype.(cdata)...}, Vector{Int}}() - indsdict!(dicts[rn+i], cdata) + cdata = tabcols(data, symbs) # Tuple(Tables.getcolumn(Tables.columns(data), x) for x in symbs) + dicts[rn + i] = Dict{Tuple{eltype.(cdata)...}, Vector{Int}}() + indsdict!(dicts[rn + i], cdata) else dicts[rn+i] = Dict(1 => collect(1:rown)) #changed to range end - sn[rn+i] = length(dicts[rn+i]) - q[rn+i] = size(rz_[i], 2) + sn[rn + i] = length(dicts[rn+i]) + q[rn + i] = size(rz_[i], 2) csp = covstrparam(repeated[i].covtype.s, q[rn+i]) - t[rn+i] = sum(csp) - tr[rn+i] = UnitRange(sum(t[1:rn+i-1]) + 1, sum(t[1:rn+i-1]) + t[rn+i]) + t[rn + i] = sum(csp) + tr[rn + i] = UnitRange(sum(t[1:rn+i-1]) + 1, sum(t[1:rn+i-1]) + t[rn+i]) updatenametype!(ct, rcnames, csp, schema[rn+i], repeated[i].covtype.s) - # emap append!(emap, fill(rn+i, t[rn+i])) end diff --git a/test/devtest.jl b/test/devtest.jl index d0650e4..d70c114 100644 --- a/test/devtest.jl +++ b/test/devtest.jl @@ -7,7 +7,7 @@ using DataFrames, CSV, StatsModels, LinearAlgebra, ForwardDiff, ForwardDiff, Opt using BenchmarkTools path = dirname(@__FILE__) cd(path) -df0 = CSV.File(path*"/csv/df0.csv"; types = [String, String, String, String, Float64, Float64]) |> DataFrame +df0 = CSV.File(path*"/csv/df0.csv"; types = [String, String, String, String, Float64, Float64, Float64]) |> DataFrame df1 = CSV.File(path*"/csv/df1.csv"; types = [String, String, String, String, Float64, Float64]) |> DataFrame ftdf = CSV.File(path*"/csv/1fptime.csv"; types = [String, String, Float64, Float64]) |> DataFrame ftdf2 = CSV.File(path*"/csv/1freparma.csv"; types = [String, String, Float64, Float64]) |> DataFrame @@ -342,3 +342,20 @@ BenchmarkTools.Trial: 2 samples with 1 evaluation. Memory estimate: 8.31 GiB, allocs estimate: 365031. =# + +#= +lmm = Metida.LMM(@formula(r2 ~ f), spatdf; +repeated = Metida.VarEffect(Metida.@covstr(x+y|1), Metida.SPEXP), +) +@benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15 + + +spatdf.ci = map(x -> CartesianIndex(x[:x], x[:y]), eachrow(spatdf)) +function Metida.edistance(mx::AbstractMatrix{<:CartesianIndex}, i::Int, j::Int) + return sqrt((mx[i, 1][1] - mx[j, 1][1])^2 + (mx[i, 1][2] - mx[j, 1][2])^2) +end +lmm = Metida.LMM(@formula(r2 ~ f), spatdf; +repeated = Metida.VarEffect(Metida.@covstr(ci|1), Metida.SPEXP; coding = Dict(:ci => Metida.RawCoding())), +) +@benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15 +=# \ No newline at end of file diff --git a/test/test.jl b/test/test.jl index 54bc17c..afcf3be 100644 --- a/test/test.jl +++ b/test/test.jl @@ -818,7 +818,6 @@ end @testset " Experimental " begin - io = IOBuffer(); lmm = Metida.LMM(@formula(r2 ~ f), spatdf; repeated = Metida.VarEffect(Metida.@covstr(x+y|1), Metida.SPEXP), @@ -827,6 +826,19 @@ end @test Metida.m2logreml(lmm) ≈ 1985.3417397854946 atol=1E-6 @test Metida.dof_satter(lmm)[1] ≈ 10.261390893063432 atol=1E-2 + + spatdf.ci = map(x -> CartesianIndex(x[:x], x[:y]), eachrow(spatdf)) + function Metida.edistance(mx::AbstractMatrix{<:CartesianIndex}, i::Int, j::Int) + return sqrt((mx[i, 1][1] - mx[j, 1][1])^2 + (mx[i, 1][2] - mx[j, 1][2])^2) + end + lmm = Metida.LMM(@formula(r2 ~ f), spatdf; + repeated = Metida.VarEffect(Metida.@covstr(ci|1), Metida.SPEXP; coding = Dict(:ci => Metida.RawCoding())), + ) + Metida.fit!(lmm) + @test Metida.m2logreml(lmm) ≈ 1985.3417397854946 atol=1E-6 + @test Metida.dof_satter(lmm)[1] ≈ 10.261390893063432 atol=1E-2 + + lmm = Metida.LMM(@formula(r2 ~ f), spatdf; repeated = Metida.VarEffect(Metida.@covstr(x+y|1), Metida.SPPOW), ) From a5d6a61c0ae77cb9f3d9237b7a21cc75b8424c09 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 26 Jan 2024 02:29:16 +0300 Subject: [PATCH 3/7] fix --- docs/src/custom.md | 2 +- src/varstruct.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/custom.md b/docs/src/custom.md index 72fb0a0..7789d09 100644 --- a/docs/src/custom.md +++ b/docs/src/custom.md @@ -163,7 +163,7 @@ Make vector of `CartesianIndex`: spatdf.ci = map(x -> CartesianIndex(x[:x], x[:y]), eachrow(spatdf)) ``` -Then use new collumn as "raw" variable with [`Metida.RawCoding`](@ref) contrast and fit the model: +Then use new column as "raw" variable with [`Metida.RawCoding`](@ref) contrast and fit the model: ```@example lmmexample lmm = Metida.LMM(@formula(r2 ~ f), spatdf; diff --git a/src/varstruct.jl b/src/varstruct.jl index cf1c3ff..49a426a 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -5,7 +5,7 @@ import StatsModels: ContrastsMatrix, AbstractContrasts, modelcols """ mutable struct RawCoding <: AbstractContrasts -Contrast for CategoricalTerm to get collumn "as it is" for model matrix. +Contrast for CategoricalTerm to get column "as it is" for model matrix. """ mutable struct RawCoding <: AbstractContrasts end From be9352f7a9e673c636fa307200415e8c7597d102 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 26 Jan 2024 02:50:19 +0300 Subject: [PATCH 4/7] try --- Project.toml | 2 ++ src/Metida.jl | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 86d8102..4b044b1 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ authors = ["Vladimir Arnautov "] version = "0.15.1" [deps] +Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -19,6 +20,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" [compat] +Compat = "3.46, 4" DiffResults = "1" Distributions = "0.21, 0.22, 0.23, 0.24, 0.25" ForwardDiff = "0.10" diff --git a/src/Metida.jl b/src/Metida.jl index 70685a3..cc1a912 100644 --- a/src/Metida.jl +++ b/src/Metida.jl @@ -3,7 +3,7 @@ __precompile__() module Metida - +using Compat using ProgressMeter, LinearAlgebra, ForwardDiff, DiffResults, Random, Optim, LineSearches, MetidaBase#, SparseArrays#, Polyester#, LoopVectorization import StatsBase, StatsModels, Distributions From 28b6467d4d3ef17aa836e408fc1b8899bbb4d7f7 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 26 Jan 2024 03:05:30 +0300 Subject: [PATCH 5/7] try --- src/varstruct.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/varstruct.jl b/src/varstruct.jl index 49a426a..8a2c574 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -263,7 +263,7 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure if isa(random[i].model, ConstantTerm) # if only ConstantTerm in the model - data_ - first is collumn (responce) data_ = data[[first(keys(data))]] else - data_ = data[StatsModels.termvars(random[i].model)] # only collumns for model + data_ = data[Tuple(StatsModels.termvars(random[i].model))] # only collumns for model end if isa(random[i].covtype.s, ZERO) schema[i] = InterceptTerm{false}() @@ -306,7 +306,7 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure if isa(repeated[i].model, ConstantTerm) # if only ConstantTerm in the model - data_ - first is collumn (responce) data_ = data[[first(keys(data))]] else - data_ = data[StatsModels.termvars(repeated[i].model)] # only collumns for model + data_ = data[Tuple(StatsModels.termvars(repeated[i].model))] # only collumns for model end schema[rn + i] = apply_schema(repeated[i].model, StatsModels.schema(data_, repeated[i].coding)) From 40d24da005ea223dd72435b3c5a0df7e0861e2f2 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 26 Jan 2024 04:38:18 +0300 Subject: [PATCH 6/7] drop 1.6 --- .github/workflows/Tier1.yml | 1 - docs/src/index.md | 4 ++++ src/Metida.jl | 1 - src/varstruct.jl | 4 ++-- 4 files changed, 6 insertions(+), 4 deletions(-) diff --git a/.github/workflows/Tier1.yml b/.github/workflows/Tier1.yml index d9481f6..7639b18 100644 --- a/.github/workflows/Tier1.yml +++ b/.github/workflows/Tier1.yml @@ -30,7 +30,6 @@ jobs: strategy: matrix: version: - - '1.6' - '1.8' - '1' os: diff --git a/docs/src/index.md b/docs/src/index.md index c38faa2..a364593 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -38,6 +38,10 @@ Implemented covariance structures: Actually Metida can fit datasets with wore than 160k observation and 40k subjects levels on PC with 64 GB RAM. This is not "hard-coded" limitation, but depends on your model and data structure. Fitting of big datasets can take a lot of time. Optimal dataset size is less than 100k observations with maximum length of block less than 400. +!!! note + + Julia v1.8 or higher required. + ## Contents ```@contents diff --git a/src/Metida.jl b/src/Metida.jl index cc1a912..a86b214 100644 --- a/src/Metida.jl +++ b/src/Metida.jl @@ -3,7 +3,6 @@ __precompile__() module Metida -using Compat using ProgressMeter, LinearAlgebra, ForwardDiff, DiffResults, Random, Optim, LineSearches, MetidaBase#, SparseArrays#, Polyester#, LoopVectorization import StatsBase, StatsModels, Distributions diff --git a/src/varstruct.jl b/src/varstruct.jl index 8a2c574..49a426a 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -263,7 +263,7 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure if isa(random[i].model, ConstantTerm) # if only ConstantTerm in the model - data_ - first is collumn (responce) data_ = data[[first(keys(data))]] else - data_ = data[Tuple(StatsModels.termvars(random[i].model))] # only collumns for model + data_ = data[StatsModels.termvars(random[i].model)] # only collumns for model end if isa(random[i].covtype.s, ZERO) schema[i] = InterceptTerm{false}() @@ -306,7 +306,7 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure if isa(repeated[i].model, ConstantTerm) # if only ConstantTerm in the model - data_ - first is collumn (responce) data_ = data[[first(keys(data))]] else - data_ = data[Tuple(StatsModels.termvars(repeated[i].model))] # only collumns for model + data_ = data[StatsModels.termvars(repeated[i].model)] # only collumns for model end schema[rn + i] = apply_schema(repeated[i].model, StatsModels.schema(data_, repeated[i].coding)) From 1ce746b2878b4caaca04ae91052ab66af02ab7b4 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 26 Jan 2024 16:36:38 +0300 Subject: [PATCH 7/7] micro opt --- Project.toml | 2 -- src/fit.jl | 4 ++-- src/gmat.jl | 11 +++++++---- src/reml.jl | 2 +- src/rmat.jl | 3 ++- src/utils.jl | 4 ++-- test/profile.pb.gz | Bin 0 -> 11360 bytes 7 files changed, 14 insertions(+), 12 deletions(-) create mode 100644 test/profile.pb.gz diff --git a/Project.toml b/Project.toml index 4b044b1..86d8102 100644 --- a/Project.toml +++ b/Project.toml @@ -6,7 +6,6 @@ authors = ["Vladimir Arnautov "] version = "0.15.1" [deps] -Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -20,7 +19,6 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" [compat] -Compat = "3.46, 4" DiffResults = "1" Distributions = "0.21, 0.22, 0.23, 0.24, 0.25" ForwardDiff = "0.10" diff --git a/src/fit.jl b/src/fit.jl index ec3b85b..a6d0788 100644 --- a/src/fit.jl +++ b/src/fit.jl @@ -63,7 +63,7 @@ Fit LMM model. * `refitinit` - true/false - if `true` - use last values for initial condition (`false` by default) * `optmethod` - Optimization method. Look at Optim.jl documentation. (Newton by default) * `singtol` - singular tolerance = 1e-8 -* `maxthreads` - maximum threads = num_cores() +* `maxthreads` - maximum threads = min(num_cores(), Threads.nthreads()) """ function fit!(lmm::LMM{T}; kwargs...) where T @@ -87,7 +87,7 @@ function fit!(lmm::LMM{T}; kwargs...) where T :refitinit ∈ kwkeys ? refitinit = kwargs[:refitinit] : refitinit = false :optmethod ∈ kwkeys ? optmethod = kwargs[:optmethod] : optmethod = :default :singtol ∈ kwkeys ? singtol = kwargs[:singtol] : singtol = 1e-8 - :maxthreads ∈ kwkeys ? maxthreads = kwargs[:maxthreads] : maxthreads = num_cores() + :maxthreads ∈ kwkeys ? maxthreads = kwargs[:maxthreads] : maxthreads = min(num_cores(), Threads.nthreads()) # If model was fitted, previous results can be used if `refitinit` == true # Before fitting clear log diff --git a/src/gmat.jl b/src/gmat.jl index c6346c9..452cf9f 100644 --- a/src/gmat.jl +++ b/src/gmat.jl @@ -148,7 +148,7 @@ function gmat!(mx, θ, ::TOEP_) if s > 1 for n = 2:s @inbounds @simd for m = 1:n-1 - mx[m, n] = de * θ[n-m+1] + mx[m, n] = de * θ[n - m + 1] end end end @@ -177,8 +177,9 @@ function gmat!(mx, θ, ::TOEPH_) end if s > 1 for n = 2:s + @inbounds mxnn = mx[n, n] @inbounds @simd for m = 1:n-1 - mx[m, n] = mx[m, m] * mx[n, n] * θ[n-m+s] + mx[m, n] = mx[m, m] * mxnn * θ[n - m + s] end end end @@ -195,8 +196,9 @@ function gmat!(mx, θ, ct::TOEPHP_) end if s > 1 && ct.p > 1 for m = 1:s - 1 + @inbounds mxmm = mx[m, m] for n = m + 1:(m + ct.p - 1 > s ? s : m + ct.p - 1) - @inbounds mx[m, n] = mx[m, m] * mx[n, n] * θ[n - m + s] + @inbounds mx[m, n] = mxmm * mx[n, n] * θ[n - m + s] end end end @@ -213,8 +215,9 @@ function gmat!(mx, θ, ::UN_) end if s > 1 for n = 2:s + @inbounds mxnn = mx[n, n] @inbounds @simd for m = 1:n - 1 - mx[m, n] = mx[m, m] * mx[n, n] * θ[s + tpnum(m, n, s)] + mx[m, n] = mx[m, m] * mxnn * θ[s + tpnum(m, n, s)] end end end diff --git a/src/reml.jl b/src/reml.jl index 5b8500c..10c29e5 100644 --- a/src/reml.jl +++ b/src/reml.jl @@ -33,7 +33,7 @@ end ################################################################################ # REML without provided β ################################################################################ -function reml_sweep_β(lmm, data, θ::Vector{T}; maxthreads::Int = 16) where T # Main optimization way - make gradient / hessian analytical / semi-analytical functions +function reml_sweep_β(lmm, data, θ::Vector{T}; maxthreads::Int = 4) where T # Main optimization way - make gradient / hessian analytical / semi-analytical functions n = length(lmm.covstr.vcovblock) N = length(lmm.data.yv) c = (N - lmm.rankx)*log(2π) diff --git a/src/rmat.jl b/src/rmat.jl index 04fcaf7..7681c7e 100644 --- a/src/rmat.jl +++ b/src/rmat.jl @@ -309,8 +309,9 @@ function unrmat(θ::AbstractVector{T}, rz) where T end if rm > 1 for m = 1:rm - 1 + @inbounds mxmm = mx[m, m] @inbounds @simd for n = m + 1:rm - mx[m, n] += mx[m, m] * mx[n, n] * θ[rm+tpnum(m, n, rm)] + mx[m, n] += mxmm * mx[n, n] * θ[rm + tpnum(m, n, rm)] end end end diff --git a/src/utils.jl b/src/utils.jl index 2bab8fa..3a8efad 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -204,8 +204,8 @@ end function logreml(lmm) -m2logreml(lmm)/2 end -function m2logreml(lmm, theta) - reml_sweep_β(lmm, LMMDataViews(lmm), theta)[1] +function m2logreml(lmm, theta; maxthreads::Int = num_cores()) + reml_sweep_β(lmm, LMMDataViews(lmm), theta; maxthreads = maxthreads)[1] end ################################################################################ diff --git a/test/profile.pb.gz b/test/profile.pb.gz new file mode 100644 index 0000000000000000000000000000000000000000..fd252a96d5b5be9e6bc12e3cd8c06cafffd88bd1 GIT binary patch literal 11360 zcmV-mET7XKiwFP!000003hjM)TohOG`0s(jh|zR3F>hnC8PRVXjiNM)$A0(lKttjc zG3HqIJ39`ojt(WS5ycBnyl;M0H8XV2pwYN~ zb|b(a)zjUt->dso)q8aNlu5eOQ}C!X-oq!>BPV5bts%Jw>+c(Mx7OgKyS47* z{Zu=C=ILsa7SE`Yk$$k|p@JK}Cub)5O+Wc)hga(D__`zZ&Z(QaeQMn+1+)AnkDXL^ z)%J;Xuk4>yH~##@x+{_<)y>~JsqX&ulj<%`nN&9|70&xV+2J*LIDm!I0E{{bpx+|^ z50(y2^qM<-NTS!{IhWSg*svu&(QDZO084HoNPL)6cHWxdiC&BD0N9;gFp$x&700#> zwf}k?eLtt{z#ZGhdo388Hs0&T3II=20n9!JAQM6Mk;g;WxE%O(_6z{(FEis_PT9K% zsF$S^@^i|bodO_d=jcSQN*?SFm8N_6Oz_xsJHvn4>|r(Y9}la!ckb{{59F_|9ldK! z?d&nBwFhqAUwzckS)U?@W^F!LvxITzk9OVJd~jIJK}$+G$HuSW;24rMdb-wZN8t#pR@ zWU;u%yO6#CDQA$rd*fCmdJUZnVA|2*L1YP$WjgcO%nqiKEZ2@-jfZvXfxe4-r@Ue2 z_erH?%qNj8&4$wp8{Z6^eM^d(l@KzElQqp^K9kw_mHXD$So~W>^5c@CSEDGQxtJ}20~N&TVm z-ss`%##W$cQW}8ctAN0r2Z{rgbxMjbOUW!Hvy{x=CYe+(8$JY*p(`taq7!cj;-&Sn z;zi~@=9BIWFlU@!OR3K}%7DhL1~MPLadBqNO2?Sb0X8fVus=9#6l^8YlY|ASEs$Sg6gzm&{UGT%5d$;I9pSWM|F zb+M&nmXcXY=6^r2ct$Wjc<`I{K93Sp%cuV~81vJ}-+C=jU0QE=V|A<;xdq=^D#DNx z^nFV;e&wpX%!Gn-<%Dm2E2s{nkJz5~)02XmOlaa;O@*J5;|P4;Q;ps}@Y=Znd-A^) z3Y>q;_d^xo#|ipYRK?v$ztV3?@~{LZ{N&+VLq#m%1bwThwr?D|IBUnu>uF5cheqX4 z9qu?n%j*NLkIue!>--*>Eb_*v?%vam(y$ z*%Rt0pRcCE@0@W2(B$#W4Tnedn|oolePh?(K1Yyvz(Tj6)qf^d)B*ZgcUlo3YjU#~WC-0a$XWygs zXNEIZs_~<2905>$^Y*Z*na8hP33Y9PAAI8opvkzLhbyx-pI`WvvPn&q!_K+mhRVey zua)2g;pXrO$G2wh&OG|9aM1%!XOFu%Y}e(n_axz#r|&l^`eX(t2sdN4KRTAQ?8NT! z!fsDEO*wvX&eA6n4}2vY@&r+(jEK3r`9bv36-;=DT_xPRaQW!9Q=21%$@uV+O(So{ z^?#Bq2=~0e8psqmK_E+^E-s?LPwjC8z_^Jq)BB%(nsWbd$|k-l{Q4e80No_4xb<+s zm>qjBy6F^sXpa*FL`P<3k4XLb@nA5_nF1&>IRs7+?ylW;Vs}bR^4V-AEU5^1Z#uSM z_{MQFuZTi|w{KGwy(oYagvPsfWIUR2>GHf;OxTPx&l}ip-Ti(CR(vkZ^>%e3^#E1; zIweN{sE(L0V&L2rR~FB6LzOn7a)QwRxxA6L_AVWL6s-EZC)}l}gcF3j(`gcT55Kvq;Z0yAy`%h&8FCUSBR>nDtcEcCkVZcj=!6Adee->?+b$&G(I?Y?zOCKvkN9M zAq^ef_{$fj#*7`XvW~FQ6EpSyMv%=m>)jsV)sjGH@f#+_yR2f4Mu zZ&Pvv*w2HE1M{O3Clt8Z4SwyDBLH0Y#U~t#$;}=R>e>x{@{%KX?odDvPUV~+(4&~L zYWcQS(f?F)f*_{MQT;B>&m7us2@}?Mf(yE}>D0ABx2N51DU@|PbLr^!Fbe(04L?+L zl$v0i^nMsl5N_uDlpG&3AoFo8VUQ0DdinIjI~jN9j{!xB^?{yfzal4ys=X!t>XDs` zMs^Wac)%1yFCydw0i(nlnenp{24DDEXbxvnhy1krWd6a&yl|(-d{)k`4QJh0OaJM{ zw>!M>nsZ0@#;l9J^_h@e2FAdzhjIik#@1W&?#CV0G?lcCw=I^;p28jy5fjmJLL!f$8qB)Zy0%N*N}~FT7+NGpD{6PJ$`3z{E}IBzY+#}!L&*#7bl3GGGR!+r75%f4PZh*(O43{9Lo^^ zW5;Jdj9a+%+SS^wG}8<4I6)YDMZv6zai`a&9AZMMH!zmAl5m2k=$odmcrfwqnPB0P zk8d3n#k-tf5wt6+rp~)D2 zWmnx)oH7_r5c*k`wQ|kelc|d_%JcM1Vk-LTA14U6Pi6I6b0>Gf4H)J$MqaTA`8fj~ z#%uy*yMTe&f-@7d_fOn8@eBJj?RI+L$+g4R3ihlB6wac-+~2Zdz?189)&~g5XfVmU z4lLa^{#IHZ1l)KYdop!I?&H{L>zObGn;*VEYy8u^4adG0@?OOEPgS%C%`ZBvEot45 zOV`piEH?>HpQSyN`E!EM_~PtwkEhO_{75bA_5kamIshjKHz)4dbac_Zg@;0fe7E_; z7kEjZTIK{{nE2c?X?cs!Zh%DKYB?B&$_kty#>VNhh7bQKe#tw+)v_R7`ZhBs2#=nh zdp>zp%-A^-!E|uCnZ07u#>LmJP8-UEaSX_#PeJkd4)PNB<(^B<8MH)=;crp$=tZ5J zAjZz~*51iHy=7Jz3~D}9(LXAjsSXJJALbfe&bfpYU!%R^j39F5IWeCxG({tfkj^lGP=0B+DQ`?*%O83IDMUT#ocmU*Gc^k%`{4{g_%0)*WR})eyA$oTYZO_ab zd>Djqr!r8FpOEJW0F#|lk_r}Y*?ryxlTTGt(clE3gZ=r-4sVP-y9{R4v1oULkEi7y z$&H=bQP_sl7L|NBL1=vI;;y7S17?iw$Alvut}6WrNR9qwgcEd~P5i(%=-*2?K?r4r z13`EfVTJz47^tn&s{8=M4rA_Cy?*tEScy|4Uaxqe>em2zczhnKIy7 z+MX?-5VOm`R7U^2!wEtcGiRP&^5jYE%r35W z{OAkFs-dF)Z{h?o#~rgFG5g7dWm{ZYwN}x`p!tP=H*=_=qG}%}2*W97UnQ@pl*kD} z7jYNvZa+L`(Du3*5NcG8{{`iaBUX;C%pH|?HmhKnOE;gX=>NGmLBK2h{=$)`E??|l zN9YWl-bkN5d0pPL6ELm!d`8ZCFJ;4_+3Ck_wii+{XU#}ISd$|F;jBy_IzMyM_RLPM zIMX&^P7t0ybm>vX^@7g&5X-*J`UH>?5 z?A}x9cUL(q+NFe4RYfITP7rQxJ0Cwc?N)S51q{zWRMCIja)KDFXWd(}V%yk#;2V<4 zfS%B6K{>&HM^Bons0Pgm!p-EH3k&j-Gh#BBP~$}w-&93!Bjg03@zG%y?+h6@^y)Dt zoKZpJTO&^mTypu&>Cc6=`0Vwe`&XUJ**hdy*yHK@or=B$!wF&~DeLIQ!DluleJsRa z)V4Whcy9XclM}#sTq+0WoAOuBPhTEAbR`oeVYD%S%*73llkZMm%7pwPo)Qv3{A>v) z2sh=Jc9=If0IG^}g6KmE%cKf1Hd0Z^o)d%~a+4lh88v9lG{~jv06$}<=Ix4)KQYo< z_={V9{Gp28-pUCAo;%am=1XcFoEf)aotrPA+C3)-trUrzqE$_mV-12kjzDeV^20M?_YBGNcb%2!6TqC{ zi;GMy2K%9E`pHM&^$9}bheP{aoKkRW>Om&NmVvwd9^sro5bj>umJ>fScf`68Ovvy6 zOQg57bArIgjO4lf2dul$?=Qksa%#}9aLcl5^ct}ce&cESmw zTZ~)%c>Kxe!DWRxSlp*|8BP!@|L2nvwhf!S_@BZJoWqaiTzsiNhnvlvLZPM=i%zFHTh+0zo)4B^=G^v zgUj*W(1%?8JjwV;rk=nwY2MI>T=e~hf200Zli*#usv=Xb!mFsgH8MAAzpD`JF2;8; z^^cm*YJaHAWh?W_uFohZ=8KG9VCq3k^Lg#}l~t`Hj2~g@pP44s2YOQkJ8$u~+u2@NQgX{?hOwp{7?$xDKz{fA9z46_by?4 zF;h=tn%iDLy22jb;qRy`X=2L&x8%ChU-@5Majl}PSTKHxsrNI@B_HTY5i|R7KXqNr zSj0%-3Fk0AhpA^XO}+;-R}_-}$^QuypZ5geDQYhou2Dy50^p(?eb?jlpyO>Qe#O+a zm+`&Gvnw7jpj^8ajxm0e^kSO_NI+32{Tu(AyACKaiWosyd5Y}Kd#R%|dy#s@vWlJhnF916 zPmm`VKY<)x%K$Y+Ub2Sq)x_bEULX>=STUaQ@pL^JbwCl~rZJwz)Z3Y6b{V*#$Wjk8 zewfsCXSv$HSC(307>{ADtbbSO#R9mfZlH-S597;4``r+=F!dg$x#XalKluw#%a0KysSkvwOjHMQJpa@mSRJC&*<*6`>}t2}W_JGPIW0(XtuOCJdvI zYjV(4$hE%!xr#Wv3-97e?)%C(GMVvYGJs1s5h|ACk>F+;UudmZu{_TBanje#=-L!( zr<`-wK|WH{UeXy)FJeoI$ldOUlh6d@#dAH^yGu!)+6(PwM;W*-2O_5#Kh4w|m}ZU# z@J6l}Luxx1XR8HxS&q1!v)&$5S49bJDdS6tmHQarB(H{r@=$l*DldjRr-O?gV2_Ht zco^fun7X_2v8sw?AT8Ir)2WE7FEf6bFsTF5QLM_>=C#$IYN`U-idgy};|H1g5Ys$F zAuIBZEXK3Qm7h^T_ljuCIUlS-tC!ce3RS5&j4Bt|TJjd45-ZWh_&QFa0#2_krSQdIwq~7xwpGWacPrzSMto7x-V87$gb0`WIKky&i9jjal zn!)%Crk=&XW~wMMueprRB?DXd0UDTGhqLhS)mBXn7(uZltK&L%$FErU{g8j?t}Aj? zID_#Fa=vM3UWzrFCcFs*y2nu#axVP`-ccH|0a)EN4w>fVPIyVrnDt1S; z;z1CCj6{`?Q{J>9@`{2}Q=XzLQ>UgF(jZuc83>@u<$&X;bC9eFdpVEAbEAN5pM(p+P#d^u40kwOa7(0 zgXVe};G4V+Y82z6ii9|Fc`b+W9AbB8jHL^AZ2g}9-raiSD(XJQ_mOt5!w5xD;ECjs z5U%b+2FVM>P04zWp_JvaK%p-gjOfS<#ZE;lz&t^*z_Fh3^)y*s$Ij)ogirVF;JeDcfTpG zAE~&?-3S$HkcCD#1RacGe}UZ&4x$~%CDe+%qPtSaJ4OnFqd^$xDmE4lXM8yM+994W zhN7_2kT+BZYL--lG2}Aq$NXb=V&vtZxB1)dhN!4E+PIA-i7JZv(gMa8ke;QY{maX# zb`=4_kCg+Y6br)bc{{M72Z*YoCf3R6g+Ed({4HmEc@dAK*fjPJ{tpNjHsKVm*p60? zmvf)B6dR)iF1X&Uq}VoVPq)7XHInN(&RpaWrsncKm0T%lpbj@i7j+!q&|r4aP9Z5+t^2 z(KPfct7Nf;Mp<g%;A8NWC;$BF!~2X8LS3d3!Tx}TxYcp(YWyJ*L8rm18i%F)Y}YZleK-6WN8&?qLYsH&Kf(aoxZz+v5h+EET9K%3|%ZbOP^pXov^9rC8j(FTRG{NjZ0v?Zpi%0k)A#XjsOX_kK@fw&CwgnrsK`ZTIw$3I>dO;yJ|N4<81ISDpVY39c3#qRao3zz|R1G>~p|GtH zpr?i|2BX2&2Ont=CUw<;+JR>>N;WAtBtVNMT3_oQ*r~fk(uIZSx=K!>Vwf905Ssvh zutS4BxFg}>HsDt#gU#PRyzrv*I?)vUo`H_PduN?=)}9Unjw+O6m}Ju#jQ)+kB21iW zgAOG)>hmHfYO&)6*wHg_X#ltc1=&3cC_`|ER#731@LE)dsE+Wfh>IdT1A}@?`bZo$ zM6z1pAc|oIYlP0G@7_Mr2GfNp3=F01S5`oLlgCT-;NPfy8!h;y2*e9avS(z(Gx_Zw zY@%5Jp0c+Ju>sGX%V+Gu+FfUn!sIRjnpyjp^xZ9HlR45F(kHx&+4yCoWa)#g8+Uc^ zGq_b0%}W60UvJsd2pweePFr?D!d1i9QEtp(B4-F%n#rtY!+_dJ^187207Nh?!C$8OmCn!$dhY5($|+9vkF z;P$fyD+gh}JtPp{>S5Hvnbv03N|N}l7Ebl{$OMjIk^w?L@jbZMO)_CP1f~p*OltrK zD}JD88!p+po5QpQ$QUdpol$K3DiKq+u6}ya2T?GXqRc%d@f`zVhfdW#zCka(9j?<` z%xbB31VqX>=lbbgbER$uTldH=p=~9bAx!7j!x(ClthUfFD<0rfeA{M@sB47W&ZxHP zA;uHGg=EW#VW=I$P!|mCcS9qsmQb2RLwjJ}6>7DG0lPz^ApZ(=gw>%)h2uAlg^YvY zXf&+^4outCDD@Wq^eW78_5&;;wz%{wzFXuH@K7v6_dfsi*`$)rE-;gYb}^fYdld~P zy)iOO(ppKs#8$57nx;){eIg`tS53#D4sHCb=B_qka;-Y*rcLb{hrjH~DgUY^wtWpa z8qGQkxDZX1Q*XpSyi(|no%(eb-+^wrYPEoZQESwp2ECf>@Rtot31p`VGdnPQySiP= z?HyMa#W5^OX9@LnyBz8ZOvSj?5B;aLuicRAs3Gt*oBIA(N8KqBq>ew@>B6Nj`mW7l))&rxtR;4M4RQh3C7Emx&s$M6Zz+ctYbnn( zzlcgjV2ly6Zq(`OjT(p%PHrhK+L>iW{uI*OIiv+RvJlAvxxJs<`nf6@Y_bJ5aQWd> z(w=Uv#xKf>J~9*>w?xj!YA4&M11&@l3Cd8nhl4kPfyf-a#Sme$IvtEPREZtA55&!3 zX@!JU@7YWAtRC1%{1=etLXVvXfsLv+0!8U-r(6l01Ukuc(BwEcuex;EQHn5jj5Ohu zj#xmj+NeTg{jnBg!m3{2ULnQlE{BX;#)TtZW>!BT`ZU%_lESO}7B3x)|Fd4#i#tz;4Sm9m^ zqt0rjbfkXU7S6VK9?-UPryy~bdwSt7QYq6>EA?-PKU$b0`m_S0fyb#v8W_?NtMx%W zx^%%PyRg@eZn>(qQAo?qkdi`16`*}aM)blgD1fpcaoWou83Oq-I4(j=iDZ{9&qvDz z`)457;1>~Wpd|`|ie@JMjYIU^CA8I+&05>>u= z_z^h2uY=pQZ2xsgXvnuAor2m1h@mD(Gr^nFTT-M2N6T>XWS@)1uk)Pgm6&mOT#OgLar4qK{t9O z{jS|A8b@}al8oF|EYjSWmAMZ|>6|+fT_IzVew3M0qzOg87(CUYBF4Wxt2@<6Z+0yB zQNs34b&C{k)LMH0P=je@Yp zRlK#tPhV6Svd8n`6wkloN(*LO6{?4eMl95aLuy%Htn)J6McNGZG`E8MEcTVio@t*h zjW9lI+EkgvHx1YIb{J%x!i>I9Kc30q!X96qyG0ns+u6R>8Emc0mS#rDVr$;VzK|Fk zCRXieGlm#}M=i{zFavr=NYM?tE=H+Eh{&IN=PmGyxJ^Ra3%V6GFLH}^k_!u6;0Co2 z`e{utLA_c%2Pa_*^yeL6j-X(*f#Zj>v=VHJf|wvmYGW|~LkbkN#fo-8E3FI~0s&FU8tVK9LuVjt~agQ4o3qoAUdwKKJX)`5y ztWX#Z3iGwDr&L5M#EMRxuo}8y#khPgXr{G`-Y$5&5ETpJ&R0(EreNTfaOzJI$ypE2$fQ#1YZG{6=rUAN` zA<(vJVH|^AvWn#u7sYaDMRp82zXlA#BaPr!4B__LioIy1DjeZ0rJokRex3&>ujn>K z7V$`mOOAuX2iAIV^&idO7#TuywhhKMELxE0pMo1~N zVtK{cE3SiJOjw?joy4nsHLqt^4O~$v>`AOdD=LJjqw6rSEQKwix784CFzPI#pp7(H zbzLP{UBcY2&@>N6@RyR#VxMWACq-w3)3qv zYUO*QHZlqbe;=Gr5Q?dV1$+-^c8EQZqc}Z05)vLQCV}ExAg&1P2({#uMb=%cVwZdb zL}U&g8%!`Qh`%dx-DI{`cZs#lDUoE|ZMzSs|-y&~1+`{2NvGZ~VN_ zy@fd{ghISXizK&@quN@cceikzElSdhpE-VLo@tja-|0&RWH=fHWXKr_rKjz#!{UKO z{Nt-g#$vO_WO}h2O!{377OSnA7+Co8*F^w6HcKC84ADsb{p&iZ;0$K$Q<=}_ZM0u? z*0yQi;tP=z(Y*~e$E+uXX)TZmgUbH-mC9?*h#q4p$4uNTEDX1vG=o_l?W2vPi;~>! z}+3@8S`qpB67L@iBBT3uI^FQRGf zZZ?LA73BqAjBQ^{;ysM`EJ~ch5+DF(^UelP1f8)h1VxrU;#-bWx+@Aw|Ba5|)mHhs z>FeUdw91D|>eLgdv5WCTsHQhZ^s$*~6}_7Hcb8My%AS!CDw3<0YD7P3p+ms|s)=PF z;OZs{G;@FN*ku%MFp1@?4!>UQ^~9$TZcA2!j^fjI+~|BjcdrIA@=vL{5Zh5l|%KaKQ%(1UhN*bPbB zK&;;=z)x>9nxWQk&BI3@m)ae;mU iDdO9-#x2|2cV$#%gtC9%IPi}V?EeG8093W-$p8Qh3r5rc literal 0 HcmV?d00001