Fri, 08 May 2020 14:46:41 -0500
proj_nonneg!
0 | 1 | ######################### |
2 | # Some utility functions | |
3 | ######################### | |
4 | ||
4
59fd17a3cea0
Add __precompile__() for what it is worth
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
5 | __precompile__() |
59fd17a3cea0
Add __precompile__() for what it is worth
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
6 | |
0 | 7 | module Util |
8 | ||
9 | ############## | |
10 | # Our exports | |
11 | ############## | |
12 | ||
13 | export map_first_slice!, | |
14 | reduce_first_slice, | |
15 | norm₂, | |
16 | γnorm₂, | |
17 | norm₂w, | |
18 | norm₂², | |
19 | norm₂w², | |
20 | norm₂₁, | |
21 | γnorm₂₁, | |
22 | dot, | |
23 | mean, | |
24 | proj_norm₂₁ball!, | |
18 | 25 | proj_nonneg!, |
0 | 26 | curry, |
8 | 27 | ⬿, |
9 | 28 | @threadsif, |
10 | 29 | @background, |
30 | @backgroundif | |
8 | 31 | |
32 | ||
33 | ########## | |
34 | # Threads | |
35 | ########## | |
36 | ||
37 | macro threadsif(threads, loop) | |
38 | return esc(:(if $threads | |
39 | Threads.@threads $loop | |
10 | 40 | else |
8 | 41 | $loop |
10 | 42 | end)) |
8 | 43 | end |
0 | 44 | |
9 | 45 | macro background(bgtask, fgtask) |
46 | return :(t = Threads.@spawn $(esc(bgtask)); | |
47 | $(esc(fgtask)); | |
48 | wait(t)) | |
49 | end | |
50 | ||
10 | 51 | macro backgroundif(threads, bgtask, fgtask) |
52 | return :(if $(esc(threads)) | |
53 | @background $(esc(bgtask)) $(esc(fgtask)) | |
54 | else | |
55 | $(esc(bgtask)) | |
56 | $(esc(fgtask)) | |
57 | end) | |
58 | end | |
59 | ||
0 | 60 | ######################## |
61 | # Functional programming | |
62 | ######################### | |
63 | ||
64 | curry = (f::Function,y...)->(z...)->f(y...,z...) | |
65 | ||
66 | ############################### | |
67 | # For working with NamedTuples | |
68 | ############################### | |
69 | ||
70 | ⬿ = merge | |
71 | ||
72 | ###### | |
73 | # map | |
74 | ###### | |
75 | ||
76 | @inline function map_first_slice!(f!, y) | |
77 | for i in CartesianIndices(size(y)[2:end]) | |
78 | @inbounds f!(@view(y[:, i])) | |
79 | end | |
80 | end | |
81 | ||
82 | @inline function map_first_slice!(x, f!, y) | |
83 | for i in CartesianIndices(size(y)[2:end]) | |
84 | @inbounds f!(@view(x[:, i]), @view(y[:, i])) | |
85 | end | |
86 | end | |
87 | ||
88 | @inline function reduce_first_slice(f, y; init=0.0) | |
89 | accum=init | |
90 | for i in CartesianIndices(size(y)[2:end]) | |
91 | @inbounds accum=f(accum, @view(y[:, i])) | |
92 | end | |
93 | return accum | |
94 | end | |
95 | ||
96 | ########################### | |
97 | # Norms and inner products | |
98 | ########################### | |
99 | ||
100 | @inline function dot(x, y) | |
101 | @assert(length(x)==length(y)) | |
102 | ||
103 | accum=0 | |
104 | for i=1:length(y) | |
105 | @inbounds accum += x[i]*y[i] | |
106 | end | |
107 | return accum | |
108 | end | |
109 | ||
110 | @inline function norm₂w²(y, w) | |
111 | #Insane memory allocs | |
112 | #return @inbounds sum(i -> y[i]*y[i]*w[i], 1:length(y)) | |
113 | accum=0 | |
114 | for i=1:length(y) | |
115 | @inbounds accum=accum+y[i]*y[i]*w[i] | |
116 | end | |
117 | return accum | |
118 | end | |
119 | ||
120 | @inline function norm₂w(y, w) | |
121 | return √(norm₂w²(y, w)) | |
122 | end | |
123 | ||
124 | @inline function norm₂²(y) | |
125 | #Insane memory allocs | |
126 | #return @inbounds sum(i -> y[i]*y[i], 1:length(y)) | |
127 | accum=0 | |
128 | for i=1:length(y) | |
129 | @inbounds accum=accum+y[i]*y[i] | |
130 | end | |
131 | return accum | |
132 | end | |
133 | ||
134 | @inline function norm₂(y) | |
135 | return √(norm₂²(y)) | |
136 | end | |
137 | ||
138 | @inline function γnorm₂(y, γ) | |
139 | hubersq = xsq -> begin | |
140 | x=√xsq | |
141 | return if x > γ | |
142 | x-γ/2 | |
143 | elseif x<-γ | |
144 | -x-γ/2 | |
145 | else | |
146 | xsq/(2γ) | |
147 | end | |
148 | end | |
149 | ||
150 | if γ==0 | |
151 | return norm₂(y) | |
152 | else | |
153 | return hubersq(norm₂²(y)) | |
154 | end | |
155 | end | |
156 | ||
157 | function norm₂₁(y) | |
158 | return reduce_first_slice((s, x) -> s+norm₂(x), y) | |
159 | end | |
160 | ||
161 | function γnorm₂₁(y,γ) | |
162 | return reduce_first_slice((s, x) -> s+γnorm₂(x, γ), y) | |
163 | end | |
164 | ||
165 | function mean(v) | |
166 | return sum(v)/prod(size(v)) | |
167 | end | |
168 | ||
169 | @inline function proj_norm₂₁ball!(y, α) | |
170 | α²=α*α | |
171 | ||
7 | 172 | if ndims(y)==3 && size(y, 1)==2 |
173 | @inbounds for i=1:size(y, 2) | |
174 | @simd for j=1:size(y, 3) | |
175 | n² = y[1,i,j]*y[1,i,j]+y[2,i,j]*y[2,i,j] | |
176 | if n²>α² | |
177 | v = α/√n² | |
178 | y[1, i, j] *= v | |
179 | y[2, i, j] *= v | |
180 | end | |
181 | end | |
182 | end | |
183 | else | |
184 | y′=reshape(y, (size(y, 1), prod(size(y)[2:end]))) | |
185 | ||
186 | @inbounds @simd for i=1:size(y′, 2)# in CartesianIndices(size(y)[2:end]) | |
187 | n² = norm₂²(@view(y′[:, i])) | |
188 | if n²>α² | |
189 | y′[:, i] .*= (α/√n²) | |
190 | end | |
0 | 191 | end |
192 | end | |
193 | end | |
194 | ||
18 | 195 | @inline function proj_nonneg!(y) |
196 | @inbounds @simd for i=1:length(y) | |
197 | if y[i] < 0 | |
198 | y[i] = 0 | |
199 | end | |
200 | end | |
201 | return y | |
202 | end | |
203 | ||
0 | 204 | end # Module |
205 |