Module:Harmonic entropy: Difference between revisions

← Older edit

@@ Line 1: / Line 1: @@
-local getArgs = require('Module:Arguments').getArgs
+local limits = require("Module:Limits")
-local u = require('Module:Utils')
+local rat = require("Module:Rational")
-local rat = require('Module:Rational')
 local p = {}
--- return measure in cents of an interval ratio, rounded to prec decimal places
+-- Compute Harmonic Shannon entropy for an interval of `c` cents
-function p.to_cents(frame)
-	local args = getArgs(frame)
-	return p._to_cents(args[1], args[2])
-end
-function p._to_cents(ratio, prec)
-	-- ratio defaults to 1
-	ratio = u.eval_num_arg(ratio, 1)
-	-- prec defaults to nil
-	prec = u.eval_num_arg(prec)
-	local result = 1200*u._log(ratio), prec
-	if prec == nil then
-		return result
-	else
-		return u._round(result, prec)
-	end
-end
--- compute harmonic Shannon entropy for an interval of `c` cents
 -- `c`, `deviation`: in cents
 -- `ratios`: an array of rational numbers
 -- `norm`: a function of rational numbers
 function p.harmonic_entropy(c, ratios, deviation, norm)
-	deviation = deviation or 17
 	norm = norm or function(ratio)
-		return math.sqrt(rat.beneditti_height(ratio))
+		return math.sqrt(rat.benedetti_height(ratio))
 	end
+	deviation = deviation or 1200 * math.log(1.01, 2)
-	local function S(x)
+	ratios = ratios
-		return math.exp(-x*x / (2*deviation*deviation)) / (deviation * math.sqrt(2*math.pi))
+		or limits.integer_limit(200, function(ratio)
+			if math.abs(rat.cents(ratio) - c) > 3 * deviation then
+				return 1 / 0
+			end
+			return norm(ratio)
+		end, 100)
+	local function gaussian(x)
+		return math.exp(-x * x / (2 * deviation * deviation)) / (deviation * math.sqrt(2 * math.pi))
 	end
-	local function Q(ratio)
-		return S(rat.cents(ratio) - c) / norm(ratio)
+	local function weighted_gaussian(ratio)
+		return gaussian(rat.cents(ratio) - c) / norm(ratio)
 	end
-	local Q_norm = 0
+	local q_norm = 0
-	for i, ratio in ipairs(ratios) do
+	for _, ratio in pairs(ratios) do
-		Q_norm = Q_norm + Q(ratio)
+		q_norm = q_norm + weighted_gaussian(ratio)
 	end
-	local function P(ratio)
+	local function probability(ratio)
-		return Q(ratio) / Q_norm
+		return weighted_gaussian(ratio) / q_norm
 	end
 	local entropy = 0
-	for i, ratio in ipairs(ratios) do
+	for _, ratio in pairs(ratios) do
-		local P_i = P(ratio)
+		local p_i = probability(ratio)
-		if P_i > 1e-5 then
+		if p_i > 1e-5 then
-			entropy = entropy - P_i * math.log(P_i) / math.log(2)
+			entropy = entropy - p_i * math.log(p_i)
 		end
 	end