-
-
Notifications
You must be signed in to change notification settings - Fork 154
/
collocation.jl
148 lines (127 loc) · 3.21 KB
/
collocation.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
abstract type CollocationKernel end
struct EpanechnikovKernel <: CollocationKernel end
struct UniformKernel <: CollocationKernel end
struct TriangularKernel <: CollocationKernel end
struct QuarticKernel <: CollocationKernel end
struct TriweightKernel <: CollocationKernel end
struct TricubeKernel <: CollocationKernel end
struct GaussianKernel <: CollocationKernel end
struct CosineKernel <: CollocationKernel end
struct LogisticKernel <: CollocationKernel end
struct SigmoidKernel <: CollocationKernel end
struct SilvermanKernel <: CollocationKernel end
function calckernel(::EpanechnikovKernel,t)
if abs(t) > 1
return 0
else
return 0.75*(1-t^2)
end
end
function calckernel(::UniformKernel,t)
if abs(t) > 1
return 0
else
return 0.5
end
end
function calckernel(::TriangularKernel,t)
if abs(t) > 1
return 0
else
return (1-abs(t))
end
end
function calckernel(::QuarticKernel,t)
if abs(t)>0
return 0
else
return (15*(1-t^2)^2)/16
end
end
function calckernel(::TriweightKernel,t)
if abs(t)>0
return 0
else
return (35*(1-t^2)^3)/32
end
end
function calckernel(::TricubeKernel,t)
if abs(t)>0
return 0
else
return (70*(1-abs(t)^3)^3)/80
end
end
function calckernel(::GaussianKernel,t)
exp(-0.5*t^2)/(sqrt(2*π))
end
function calckernel(::CosineKernel,t)
if abs(t)>0
return 0
else
return (π*cos(π*t/2))/4
end
end
function calckernel(::LogisticKernel,t)
1/(exp(t)+2+exp(-t))
end
function calckernel(::SigmoidKernel,t)
2/(π*(exp(t)+exp(-t)))
end
function calckernel(::SilvermanKernel,t)
sin(abs(t)/2+π/4)*0.5*exp(-abs(t)/sqrt(2))
end
function construct_t1(t,tpoints)
hcat(ones(eltype(tpoints),length(tpoints)),tpoints.-t)
end
function construct_t2(t,tpoints)
hcat(ones(eltype(tpoints),length(tpoints)),tpoints.-t,(tpoints.-t).^2)
end
function construct_w(t,tpoints,h,kernel)
W = @. calckernel((kernel,),(tpoints-t)/h)/h
Diagonal(W)
end
"""
```julia
u′,u = collocate_data(data,tpoints,kernel=SigmoidKernel())
```
Computes a non-parametrically smoothed estimate of `u'` and `u`
given the `data`, where each column is a snapshot of the timeseries at
`tpoints[i]`.
For kernels, the following exist:
- EpanechnikovKernel
- UniformKernel
- TriangularKernel
- QuarticKernel
- TriweightKernel
- TricubeKernel
- GaussianKernel
- CosineKernel
- LogisticKernel
- SigmoidKernel
- SilvermanKernel
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2631937/
"""
function collocate_data(data,tpoints,kernel=TriangularKernel())
_one = oneunit(first(data))
_zero = zero(first(data))
e1 = [_one;_zero]
e2 = [_zero;_one;_zero]
n = length(tpoints)
h = (n^(-1/5))*(n^(-3/35))*((log(n))^(-1/16))
Wd = similar(data, n, size(data,1))
WT1 = similar(data, n, 2)
WT2 = similar(data, n, 3)
x = map(tpoints) do _t
T1 = construct_t1(_t,tpoints)
T2 = construct_t2(_t,tpoints)
W = construct_w(_t,tpoints,h,kernel)
mul!(Wd,W,data')
mul!(WT1,W,T1)
mul!(WT2,W,T2)
(e2'*((T2'*WT2)\T2'))*Wd,(e1'*((T1'*WT1)\T1'))*Wd
end
estimated_derivative = reduce(hcat,transpose.(first.(x)))
estimated_solution = reduce(hcat,transpose.(last.(x)))
estimated_derivative,estimated_solution
end