Skip to content

Commit 34176af

Browse files
andrewjradcliffeandrewjradcliffe
committed
Organize, minor cleaning
1 parent 150f107 commit 34176af

File tree

2 files changed

+72
-40
lines changed

2 files changed

+72
-40
lines changed

‎src/vnorm.jl‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,10 @@
77
_denom(A, dims) = prod(d -> size(A, d), dims)
88
_denom(A, ::Colon) = length(A)
99

10+
@inline _reducedsize(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M} = ntuple(d -> d ∈ dims ? 1 : size(A, d), Val(N))
11+
1012
# p-norms
11-
_norm0(A, dims::NTuple{M, Int}) where {M} =
12-
fill!(similar(A, ntuple(d -> d ∈ dims ? 1 : size(A, d), Val(M))), _denom(A, dims))
13+
_norm0(A, dims::NTuple{M, Int}) where {M} = fill!(similar(A, _reducedsize(A, dims)), _denom(A, dims))
1314
_norm0(A, ::Colon) = float(length(A))
1415

1516
_vnorm(A, p::Rational{T}, dims) where {T} = _vnorm(A, float(p), dims)

‎src/vstats.jl‎

Lines changed: 69 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -4,25 +4,10 @@
44
#
55
#
66
############################################################################################
7-
87
# Statistical things
9-
function vmse(x, y; dims=:)
10-
c = 1 / _denom(x, dims)
11-
vmapreducethen((xᵢ, yᵢ) -> abs2(xᵢ - yᵢ) , +, z -> c * z, x, y, dims=dims)
12-
end
13-
function vtmse(x, y; dims=:)
14-
c = 1 / _denom(x, dims)
15-
vtmapreducethen((xᵢ, yᵢ) -> abs2(xᵢ - yᵢ) , +, z -> c * z, x, y, dims=dims)
16-
end
17-
function vrmse(x, y; dims=:)
18-
c = 1 / _denom(x, dims)
19-
vmapreducethen((xᵢ, yᵢ) -> abs2(xᵢ - yᵢ) , +, z -> √(c * z), x, y, dims=dims)
20-
end
21-
function vtrmse(x, y; dims=:)
22-
c = 1 / _denom(x, dims)
23-
vtmapreducethen((xᵢ, yᵢ) -> abs2(xᵢ - yᵢ) , +, z -> √(c * z), x, y, dims=dims)
24-
end
258

9+
################
10+
# Means
2611
function vmean(f, A; dims=:)
2712
c = 1 / _denom(A, dims)
2813
vmapreducethen(f, +, x -> c * x, A, dims=dims)
@@ -57,6 +42,9 @@ end
5742
_xlogx(x::T) where {T} = ifelse(iszero(x), zero(T), x * log(x))
5843
_xlogy(x::T, y::T) where {T} = ifelse(iszero(x) & !isnan(y), zero(T), x * log(y))
5944

45+
################
46+
# Entropies
47+
6048
ventropy(A; dims=:) = vmapreducethen(_xlogx, +, -, A, dims=dims)
6149
ventropy(A, b::Real; dims=:) = (c = -1 / log(b); vmapreducethen(_xlogx, +, x -> x * c, A, dims=dims ))
6250

@@ -65,7 +53,7 @@ vcrossentropy(p, q, b::Real; dims=:) = (c = -1 / log(b); vmapreducethen(_xlogy,
6553

6654
# max-, Shannon, collision and min- entropy assume that p ∈ ℝⁿ, pᵢ ≥ 0, ∑pᵢ=1
6755
_vmaxentropy(p, dims::NTuple{M, Int}) where {M} =
68-
fill!(similar(p, ntuple(d -> d ∈ dims ? 1 : size(p, d), Val(M))), log(_denom(p, dims)))
56+
fill!(similar(p, _reducedsize(p, dims)), log(_denom(p, dims)))
6957
_vmaxentropy(p, ::Colon) = log(length(p))
7058
vmaxentropy(p; dims=:) = _vmaxentropy(p, dims)
7159
vshannonentropy(p; dims=:) = vmapreducethen(_xlogx, +, -, p, dims=dims)
@@ -78,6 +66,7 @@ _vrenyientropy(p, α::T, dims) where {T<:AbstractFloat} =
7866
(c = one(T) / (one(T) - α); vmapreducethen(x -> exp(α * log(x)), +, x -> c * log(x), p, dims=dims))
7967
_vrenyientropy(p, α::Rational{T}, dims) where {T} = _vrenyientropy(p, float(α), dims)
8068
function vrenyientropy(p, α::Real; dims=:)
69+
α < 0 && throw(ArgumentError("Order of Rényi entropy not legal, $(α) < 0."))
8170
if α ≈ 0
8271
vmaxentropy(p, dims=dims)
8372
elseif α ≈ 1
@@ -93,23 +82,41 @@ end
9382
# Loosened restrictions: p ∈ ℝⁿ, pᵢ ≥ 0, ∑pᵢ > 1; that is, if one normalized p, a valid
9483
# probability vector would be produced. Thus, H(x, α) = (α/(1-α)) * (1/α * log∑xᵢ^α - log∑xᵢ)
9584
# H(x, α) = (α / (1 - α)) * ((1/α) * log(sum(z -> z^α, x)) - log(sum(x)))
96-
vrenyientropynorm(p, α::Real; dims=:) =
97-
vrenyientropy(p, α, dims=dims) .- (α/(1-α)) .* log.(vnorm(p, 1, dims=dims))
98-
99-
vrenyientropy(x2n, 1.5)
100-
renyientropy(x2n, 1.5)
101-
102-
den = sum(abs2, x2)
103-
sum(abs2.(x2)./ den)
104-
sum(abs2, x2 ./ √den)
105-
106-
√(abs2(1 / sum(abs, x2)) * sum(abs2, x2))
107-
(1 / sum(abs, x2)) * norm(x2)
108-
norm(x2n)
109-
norm(x2) / norm(x2, 1)
110-
log(norm(x2))
111-
log(norm(x2)) - log(norm(x2, 1))
85+
# vrenyientropynorm(p, α::Real; dims=:) =
86+
# vrenyientropy(p, α, dims=dims) .- (α/(1-α)) .* log.(vnorm(p, 1, dims=dims))
87+
88+
####
89+
# threaded versions
90+
_vtmaxentropy(p, dims::NTuple{M, Int}) where {M} =
91+
fill!(similar(p, _reducedsize(p, dims)), log(_denom(p, dims)))
92+
_vtmaxentropy(p, ::Colon) = log(length(p))
93+
vtmaxentropy(p; dims=:) = _vtmaxentropy(p, dims)
94+
vtshannonentropy(p; dims=:) = vtmapreducethen(_xlogx, +, -, p, dims=dims)
95+
vtcollisionentropy(p; dims=:) = vtmapreducethen(abs2, +, x -> -log(x), p, dims=dims)
96+
vtminentropy(p; dims=:) = vtmapreducethen(identity, max, x -> -log(x), p, dims=dims)
97+
98+
_vtrenyientropy(p, α::T, dims) where {T<:Integer} =
99+
(c = one(T) / (one(T) - α); vtmapreducethen(x -> x^α, +, x -> c * log(x), p, dims=dims))
100+
_vtrenyientropy(p, α::T, dims) where {T<:AbstractFloat} =
101+
(c = one(T) / (one(T) - α); vtmapreducethen(x -> exp(α * log(x)), +, x -> c * log(x), p, dims=dims))
102+
_vtrenyientropy(p, α::Rational{T}, dims) where {T} = _vtrenyientropy(p, float(α), dims)
103+
function vtrenyientropy(p, α::Real; dims=:)
104+
α < 0 && throw(ArgumentError("Order of Rényi entropy not legal, $(α) < 0."))
105+
if α ≈ 0
106+
vtmaxentropy(p, dims=dims)
107+
elseif α ≈ 1
108+
vtshannonentropy(p, dims=dims)
109+
elseif α ≈ 2
110+
vtcollisionentropy(p, dims=dims)
111+
elseif isinf(α)
112+
vtminentropy(p, dims=dims)
113+
else
114+
_vtrenyientropy(p, α, dims)
115+
end
116+
end
112117

118+
################
119+
# Divergences
113120

114121
# # StatsBase handling of pᵢ = qᵢ = 0
115122
# _xlogxdy(x::T, y::T) where {T} = _xlogy(x, ifelse(iszero(x) & iszero(y), zero(T), x / y))
@@ -118,9 +125,15 @@ log(norm(x2)) - log(norm(x2, 1))
118125
_klterm(x::T, y::T) where {T} = _xlogy(x, x) - _xlogy(x, y)
119126
vkldivergence(p, q; dims=:) = vvmapreduce(_klterm, +, p, q, dims=dims)
120127
vkldivergence(p, q, b::Real; dims=:) = (c = 1 / log(b); vmapreducethen(_klterm, +, x -> x * c, p, q, dims=dims))
128+
vtkldivergence(p, q; dims=:) = vtmapreduce(_klterm, +, p, q, dims=dims)
129+
vtkldivergence(p, q, b::Real; dims=:) = (c = 1 / log(b); vtmapreducethen(_klterm, +, x -> x * c, p, q, dims=dims))
121130

131+
# generalized KL divergence sum(a*log(a/b)-a+b)
132+
vgkldiv(x, y; dims=:) = vvmapreduce((xᵢ, yᵢ) -> xᵢ * (log(xᵢ) - log(yᵢ)) - xᵢ + yᵢ, +, x, y, dims=dims)
133+
vtgkldiv(x, y; dims=:) = vtmapreduce((xᵢ, yᵢ) -> xᵢ * (log(xᵢ) - log(yᵢ)) - xᵢ + yᵢ, +, x, y, dims=dims)
122134

123-
135+
################
136+
# Deviations
124137

125138
vcounteq(x, y; dims=:) = vvmapreduce(==, +, x, y, dims=dims)
126139
vtcounteq(x, y; dims=:) = vtmapreduce(==, +, x, y, dims=dims)
@@ -139,7 +152,25 @@ end
139152
vmaxad(x, y; dims=:) = vvmapreduce((xᵢ, yᵢ) -> abs(xᵢ - yᵢ) , max, x, y, dims=dims)
140153
vtmaxad(x, y; dims=:) = vtmapreduce((xᵢ, yᵢ) -> abs(xᵢ - yᵢ) , max, x, y, dims=dims)
141154

155+
function vmse(x, y; dims=:)
156+
c = 1 / _denom(x, dims)
157+
vmapreducethen((xᵢ, yᵢ) -> abs2(xᵢ - yᵢ) , +, z -> c * z, x, y, dims=dims)
158+
end
159+
function vtmse(x, y; dims=:)
160+
c = 1 / _denom(x, dims)
161+
vtmapreducethen((xᵢ, yᵢ) -> abs2(xᵢ - yᵢ) , +, z -> c * z, x, y, dims=dims)
162+
end
163+
function vrmse(x, y; dims=:)
164+
c = 1 / _denom(x, dims)
165+
vmapreducethen((xᵢ, yᵢ) -> abs2(xᵢ - yᵢ) , +, z -> √(c * z), x, y, dims=dims)
166+
end
167+
function vtrmse(x, y; dims=:)
168+
c = 1 / _denom(x, dims)
169+
vtmapreducethen((xᵢ, yᵢ) -> abs2(xᵢ - yᵢ) , +, z -> √(c * z), x, y, dims=dims)
170+
end
142171

143-
# generalized KL divergence sum(a*log(a/b)-a+b)
144-
vgkldiv(x, y; dims=:) = vvmapreduce((xᵢ, yᵢ) -> xᵢ * (log(xᵢ) - log(yᵢ)) - xᵢ + yᵢ, +, x, y, dims=dims)
145-
vtgkldiv(x, y; dims=:) = vtmapreduce((xᵢ, yᵢ) -> xᵢ * (log(xᵢ) - log(yᵢ)) - xᵢ + yᵢ, +, x, y, dims=dims)
172+
# Match names with StatsBase
173+
const vmsd = vmse
174+
const vtmse = vtmsd
175+
const vrmsd = vrmse
176+
const vtrmsd = vtrmse

0 commit comments

Comments
 (0)