Skip to content

Commit 67901aa

Browse files
luyouepiusfluyouepiusf
committed
initial commit
1 parent 7e1452d commit 67901aa

File tree

5 files changed

+411
-3
lines changed

5 files changed

+411
-3
lines changed

‎Project.toml‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,13 @@
11
name = "IntervalCensoredMultistate"
22
uuid = "3506c8fb-fa20-493b-899a-e67d94b3cb46"
33
authors = ["Lu You"]
4-
version = "1.0.0-DEV"
4+
version = "0.1.0"
5+
6+
[deps]
7+
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
58

69
[compat]
7-
julia = "1.10.10"
10+
julia = "1.4"
811

912
[extras]
1013
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

‎src/IntervalCensoredMultistate.jl‎

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,15 @@
11
module IntervalCensoredMultistate
22

3-
# Write your package code here.
3+
import LinearAlgebra
4+
import DataFrames
5+
6+
export greet_your_package_name
7+
include("infer_transitions.jl")
8+
9+
export reg_est_bone
10+
include("reg_est_bone.jl")
11+
12+
export simplify
13+
include("simplify.jl")
414

515
end

‎src/infer_transitions.jl‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
function infer_transitions(possible_transition)
2+
nstates=size(possible_transition,1)
3+
4+
possible_transition_tuple = Vector{Tuple{Int64,Int64}}(undef,0)
5+
for idx_begin = 1:nstates
6+
for idx_end = 1:nstates
7+
if possible_transition[idx_begin,idx_end]
8+
push!(possible_transition_tuple,CartesianIndex(idx_begin,idx_end))
9+
end
10+
end
11+
end
12+
13+
reachable_within=fill(Inf,nstates,nstates)
14+
for ii in nstates:(-1):0
15+
for idx in findall((possible_transition^ii).>=1)
16+
reachable_within[idx]=ii
17+
end
18+
end
19+
possible_sub=.!(isinf.(reachable_within))
20+
21+
F_next_states=Vector{Vector{Int64}}(undef, nstates)
22+
F_sub_states=Matrix{Vector{Int64}}(undef, nstates, nstates)
23+
for idx_begin in 1:nstates
24+
idx_sub=possible_transition[idx_begin,:]
25+
if any(idx_sub) F_next_states[idx_begin]=findall(idx_sub) end
26+
end
27+
for idx_begin in 1:nstates
28+
for idx_end in 1:nstates
29+
idx_sub=possible_sub[idx_begin,:].&possible_sub[:,idx_end]
30+
if any(idx_sub) F_sub_states[idx_begin,idx_end]=findall(idx_sub) end
31+
end
32+
end
33+
34+
absorbing_states=.!any(possible_transition, dims=2)
35+
36+
cycles=fill(false,nstates)
37+
for ii in nstates:(-1):1
38+
for idx in findall(LinearAlgebra.diag(possible_transition^ii).>=1)
39+
cycles[idx]=true
40+
end
41+
end
42+
43+
return (
44+
possible_transition_tuple,
45+
reachable_within,
46+
possible_sub,
47+
F_next_states,
48+
F_sub_states,
49+
absorbing_states,
50+
cycles)
51+
end

‎src/reg_est_bone.jl‎

Lines changed: 252 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,252 @@
1+
function reg_est_bone(
2+
TTij,
3+
ssij,
4+
nobss,
5+
zzi,
6+
F_zidx,
7+
possible_transition,
8+
ntimepoints;
9+
niter=1000,tol=1e-4,factor=0.8)
10+
11+
nind=size(TTij,1)
12+
ndimz=size(zzi,2)
13+
nstates=size(possible_transition,1)
14+
eye_nstates=Matrix{Float64}(LinearAlgebra.I,nstates,nstates)
15+
(
16+
possible_transition_tuple,
17+
reachable_within,
18+
possible_sub,
19+
F_next_states,
20+
F_sub_states,
21+
absorbing_states,
22+
cycles)=infer_transitions(possible_transition)
23+
24+
F_nzidx=fill(0,nstates,nstates)
25+
for (idx_begin,idx_end) in possible_transition_tuple
26+
F_nzidx[idx_begin,idx_end]=length(F_zidx[idx_begin,idx_end])
27+
end
28+
29+
F_betaz_est = Matrix{Vector}(undef, nstates, nstates)
30+
F_betaz_old = Matrix{Vector}(undef, nstates, nstates)
31+
F_eps = fill(0.0, nstates, nstates)
32+
F_hazard0 = Matrix{Vector}(undef, nstates, nstates)
33+
34+
for idx_begin = 1:nstates
35+
for idx_end = 1:nstates
36+
if possible_transition[idx_begin,idx_end]
37+
F_betaz_est[idx_begin,idx_end] = fill(0.0, F_nzidx[idx_begin,idx_end])
38+
F_betaz_old[idx_begin,idx_end] = fill(0.0, F_nzidx[idx_begin,idx_end])
39+
F_hazard0[idx_begin,idx_end] = fill(0.002, ntimepoints)
40+
end
41+
end
42+
end
43+
F_SS_betaz = Matrix{Vector}(undef, nstates, nstates)
44+
F_II_betaz = Matrix{Matrix}(undef, nstates, nstates)
45+
46+
F_nonzero_g=Matrix{Vector{Int64}}(undef,nind,size(ssij,2))
47+
F_nonzero_gd=Array{Vector{Int64},3}(undef,nind,size(ssij,2),nstates)
48+
for ii in 1:nind
49+
for jj in 1:(nobss[ii]-1)
50+
a_sub_states=Vector{Int64}(undef,0)
51+
for idx_begin in findall(ssij[ii,jj,:])
52+
for idx_end in findall(ssij[ii,jj+1,:])
53+
if isassigned(F_sub_states,idx_begin,idx_end) append!(a_sub_states,F_sub_states[idx_begin,idx_end]) end
54+
end
55+
end
56+
unique_sub_states=sort(unique(a_sub_states))
57+
F_nonzero_g[ii,jj]=unique_sub_states
58+
a_next_states=Vector{Int64}(undef,0)
59+
for idx_begin in unique_sub_states
60+
if isassigned(F_next_states,idx_begin)
61+
F_nonzero_gd[ii,jj,idx_begin]=F_next_states[idx_begin]
62+
else
63+
F_nonzero_gd[ii,jj,idx_begin]=[]
64+
end
65+
end
66+
end
67+
end
68+
69+
F_one_g=fill(false,nind,size(ssij,2),nstates)
70+
for ii in 1:nind
71+
for jj in 1:(nobss[ii]-1)
72+
if sum(ssij[ii,jj,:])==1 && all(ssij[ii,jj,:].==ssij[ii,jj+1,:])
73+
ss=findfirst(ssij[ii,jj,:])
74+
if !cycles[ss]
75+
F_one_g[ii,jj,ss]=true
76+
end
77+
end
78+
end
79+
end
80+
81+
## define blocks
82+
nblock=fill(0,nind)
83+
block2jj=fill(-1,nind,size(TTij)[2],2)
84+
block2tt=fill(-1,nind,size(TTij)[2],2)
85+
for ii in 1:nind
86+
temp_nblock=0
87+
temp_block2jj=fill(-1,size(TTij)[2],2)
88+
temp_block2tt=fill(-1,size(TTij)[2],2)
89+
temp_begin=nobss[ii]
90+
temp_end=nobss[ii]
91+
for jj in (nobss[ii]-1):(-1):1
92+
if sum(ssij[ii,jj,:])==1
93+
temp_begin=jj
94+
temp_nblock=temp_nblock+1
95+
temp_block2jj[temp_nblock,1]=temp_begin
96+
temp_block2jj[temp_nblock,2]=temp_end
97+
temp_block2tt[temp_nblock,1]=TTij[ii,temp_begin]
98+
temp_block2tt[temp_nblock,2]=TTij[ii,temp_end]
99+
temp_end=jj
100+
end
101+
end
102+
nblock[ii]=temp_nblock
103+
block2jj[ii,1:temp_nblock,:]=temp_block2jj[temp_nblock:(-1):1,:]
104+
block2tt[ii,1:temp_nblock,:]=temp_block2tt[temp_nblock:(-1):1,:]
105+
end
106+
107+
prob_by_jj=fill(NaN,nstates,nstates,size(TTij)[2])
108+
prob_lt_tt=fill(NaN,nstates,nstates,ntimepoints)
109+
prob_gt_tt=fill(NaN,nstates,nstates,ntimepoints)
110+
prob_ge_tt=fill(NaN,nstates,nstates,ntimepoints)
111+
112+
F_Egd = fill(0.0, nind, nstates, nstates, ntimepoints)
113+
F_Eg = fill(0.0, nind, nstates, ntimepoints)
114+
F_expr = fill(NaN, nind, nstates, nstates)
115+
F_zexpr = Array{Vector{Float64},3}(undef, nind, nstates, nstates)
116+
F_zzexpr = Array{Matrix{Float64},3}(undef, nind, nstates, nstates)
117+
118+
start_time=time()
119+
for iter in 1:niter
120+
println(iter," - ")
121+
current_time=time()
122+
println(current_time-start_time)
123+
for ii in 1:nind
124+
prob=zeros(nstates,nstates,ntimepoints)
125+
expr=zeros(nstates,nstates)
126+
for (idx_begin,idx_end) in possible_transition_tuple
127+
zz_temp=zzi[ii,F_zidx[idx_begin,idx_end]]
128+
expr[idx_begin,idx_end]=exp(LinearAlgebra.dot(F_betaz_est[idx_begin,idx_end],zz_temp))
129+
F_expr[ii,idx_begin,idx_end]=expr[idx_begin,idx_end]
130+
F_zexpr[ii,idx_begin,idx_end] = F_expr[ii,idx_begin,idx_end] * zz_temp
131+
F_zzexpr[ii,idx_begin,idx_end] = F_zexpr[ii,idx_begin,idx_end] * transpose(zz_temp)
132+
prob[idx_begin,idx_end,:]=1.0.-exp.(-F_hazard0[idx_begin,idx_end]*expr[idx_begin,idx_end])
133+
end
134+
for idx_begin in 1:nstates
135+
prob[idx_begin,idx_begin,:]=1.0.-sum(prob[idx_begin,:,:],dims=1)
136+
end
137+
prob_by_block=fill(NaN,size(ssij)[2])
138+
prob_by_jj=fill(NaN,nstates,nstates,size(ssij)[2])
139+
prob_lt_tt=fill(NaN,nstates,nstates,ntimepoints) # can be defined outside the loop
140+
prob_gt_tt=fill(NaN,nstates,nstates,ntimepoints) # can be defined outside the loop
141+
prob_ge_tt=fill(NaN,nstates,nstates,ntimepoints) # can be defined outside the loop
142+
143+
for jj in 1:(nobss[ii]-1)
144+
prob_temp=copy(eye_nstates)
145+
for tt in (TTij[ii,jj]+1):TTij[ii,jj+1]
146+
prob_lt_tt[:,:,tt]=prob_temp
147+
prob_temp=prob_temp*prob[:,:,tt]
148+
end
149+
prob_by_jj[:,:,jj]=prob_temp
150+
prob_temp=copy(eye_nstates)
151+
for tt in (TTij[ii,jj+1]):(-1):(TTij[ii,jj]+1)
152+
prob_gt_tt[:,:,tt]=prob_temp
153+
prob_temp=prob[:,:,tt]*prob_temp
154+
prob_ge_tt[:,:,tt]=prob_temp
155+
end
156+
end
157+
F_lt_jj = Vector{Matrix{Float64}}(undef,size(ssij)[2])
158+
F_gt_jj = Vector{Matrix{Float64}}(undef,size(ssij)[2])
159+
for bb in 1:nblock[ii]
160+
block_temp=fill(1.0,1,sum(ssij[ii,block2jj[ii,bb,1],:]))
161+
for jj in block2jj[ii,bb,1]:(block2jj[ii,bb,2]-1)
162+
F_lt_jj[jj]=block_temp
163+
block_temp=block_temp*prob_by_jj[ssij[ii,jj,:],ssij[ii,jj+1,:],jj]
164+
end
165+
prob_by_block[bb]=sum(block_temp)
166+
block_temp=fill(1.0,sum(ssij[ii,block2jj[ii,bb,2],:]),1)
167+
for jj in (block2jj[ii,bb,2]-1):(-1):block2jj[ii,bb,1]
168+
F_gt_jj[jj]=block_temp
169+
block_temp=prob_by_jj[ssij[ii,jj,:],ssij[ii,jj+1,:],jj]*block_temp
170+
end
171+
end
172+
for bb in 1:nblock[ii]
173+
for jj in block2jj[ii,bb,1]:(block2jj[ii,bb,2]-1)
174+
for tt in (TTij[ii,jj]+1):(TTij[ii,jj+1])
175+
for a_sub_idx_begin in F_nonzero_g[ii,jj]
176+
if F_one_g[ii,jj,a_sub_idx_begin]
177+
g_temp=1.0
178+
else
179+
g_temp=sum(
180+
F_lt_jj[jj]*
181+
prob_lt_tt[ssij[ii,jj,:],a_sub_idx_begin:a_sub_idx_begin,tt]*
182+
prob_ge_tt[a_sub_idx_begin:a_sub_idx_begin,ssij[ii,jj+1,:],tt]*
183+
F_gt_jj[jj])/prob_by_block[bb]
184+
end
185+
F_Eg[ii,a_sub_idx_begin,tt]=g_temp
186+
for a_sub_idx_end in F_nonzero_gd[ii,jj,a_sub_idx_begin]
187+
gd_temp=sum(
188+
F_lt_jj[jj]*
189+
prob_lt_tt[ssij[ii,jj,:],a_sub_idx_begin:a_sub_idx_begin,tt]*
190+
prob[a_sub_idx_begin:a_sub_idx_begin,a_sub_idx_end:a_sub_idx_end,tt]*
191+
prob_gt_tt[a_sub_idx_end:a_sub_idx_end,ssij[ii,jj+1,:],tt]*
192+
F_gt_jj[jj])/prob_by_block[bb]
193+
F_Egd[ii,a_sub_idx_begin,a_sub_idx_end,tt]=gd_temp
194+
end
195+
end
196+
end
197+
end
198+
end
199+
end
200+
for (idx_begin,idx_end) in possible_transition_tuple
201+
AA0=fill(0.0,ntimepoints)
202+
AA1=fill(0.0,F_nzidx[idx_begin,idx_end],ntimepoints)
203+
AA2=fill(0.0,F_nzidx[idx_begin,idx_end],F_nzidx[idx_begin,idx_end],ntimepoints)
204+
BB1=fill(0.0,F_nzidx[idx_begin,idx_end],ntimepoints)
205+
BB2=fill(0.0,F_nzidx[idx_begin,idx_end],F_nzidx[idx_begin,idx_end],ntimepoints)
206+
for tt in 1:ntimepoints
207+
for ii in 1:nind
208+
AA0[tt]=AA0[tt]+(F_Eg[ii,idx_begin,tt]-F_Egd[ii,idx_begin,idx_end,tt]/2.0)*F_expr[ii,idx_begin,idx_end]
209+
AA1[:,tt]=AA1[:,tt]+(F_Eg[ii,idx_begin,tt]-F_Egd[ii,idx_begin,idx_end,tt]/2.0)*F_zexpr[ii,idx_begin,idx_end]
210+
AA2[:,:,tt]=AA2[:,:,tt]+(F_Eg[ii,idx_begin,tt]-F_Egd[ii,idx_begin,idx_end,tt]/2.0)*F_zzexpr[ii,idx_begin,idx_end]
211+
end
212+
if AA0[tt]==0.0
213+
BB1[:,tt].=0.0
214+
BB2[:,:,tt].=0.0
215+
else
216+
BB1[:,tt]=AA1[:,tt]/AA0[tt]
217+
BB2[:,:,tt]=AA2[:,:,tt]/AA0[tt]
218+
end
219+
end
220+
SS_betaz=fill(0.0,F_nzidx[idx_begin,idx_end]) # can be defined outside the loop
221+
II_betaz=fill(0.0,F_nzidx[idx_begin,idx_end],F_nzidx[idx_begin,idx_end]) # can be defined outside the loop
222+
for tt in 1:ntimepoints
223+
for ii in 1:nind
224+
SS_betaz=SS_betaz+F_Egd[ii,idx_begin,idx_end,tt]*zzi[ii,F_zidx[idx_begin,idx_end]]-F_Egd[ii,idx_begin,idx_end,tt]*BB1[:,tt]
225+
II_betaz=II_betaz+F_Egd[ii,idx_begin,idx_end,tt]*(BB1[:,tt]*transpose(BB1[:,tt])-BB2[:,:,tt])
226+
end
227+
end
228+
F_betaz_old[idx_begin,idx_end]=deepcopy(F_betaz_est[idx_begin,idx_end])
229+
F_betaz_est[idx_begin,idx_end]=F_betaz_est[idx_begin,idx_end]-factor*(II_betaz\SS_betaz)
230+
F_eps[idx_begin,idx_end]=sum(abs.(F_betaz_est[idx_begin,idx_end]-F_betaz_old[idx_begin,idx_end]))
231+
F_II_betaz[idx_begin,idx_end]=II_betaz
232+
for tt = 1:ntimepoints
233+
sum_Egd = sum(F_Egd[:,idx_begin,idx_end,tt]) # can be defined outside the loop
234+
if AA0[tt]==0.0
235+
F_hazard0[idx_begin,idx_end][tt]=0.0
236+
else
237+
F_hazard0[idx_begin,idx_end][tt]=sum_Egd/AA0[tt]
238+
end
239+
end
240+
end
241+
max_eps=sum(F_eps)
242+
if max_eps<tol break end
243+
end
244+
245+
return(
246+
F_betaz_est,
247+
F_II_betaz,
248+
F_hazard0,
249+
F_Eg,
250+
F_eps)
251+
end
252+

0 commit comments

Comments
 (0)