Hi all,
The so-called Hungarian Algorithm (Hungarian algorithm - Wikipedia) can be useful, for example, for matching instruments (given instrument range / ease of playing note within the range) to a collection of notes to play.
I couldn’t find any implementation in SC, although it’s pretty common in other languages. Am I missing something?
Thanks!
1 Like
Here is a simple implementation for reference.
(
// Subtract the minimum value from each row
~rowMin = { |m|
m.collect { |r|
var mn = r.minItem;
r.collect { |v| v - mn }
}
};
// Subtract the minimum value from each column
~colMin = { |m|
var sz = m.size, r = m.deepCopy;
sz.do { |j|
var mn = inf;
sz.do { |i| mn = min(mn, r[i][j]) };
if(mn > 0) { sz.do { |i| r[i][j] = r[i][j] - mn } };
};
r
};
// Find zeros and "star" the first uncovered zero in each row (using a new array)
~starZero = { |m|
var sz = m.size, mask = Array.fill(sz, { Array.fill(sz, 0) });
var rC = Array.fill(sz, false), cC = Array.fill(sz, false);
sz.do { |i|
sz.do { |j|
if((m[i][j] == 0) && rC[i].not && cC[j].not) {
mask[i][j] = 1; rC[i] = true; cC[j] = true;
}
}
};
mask
};
// Identify which columns contain a starred zero.
~coverStar = { |mask|
var sz = mask.size, cC = Array.fill(sz, false);
sz.do { |j| sz.do { |i| if(mask[i][j] == 1) { cC[j] = true } } };
cC
};
// Locate an uncovered zero
~findZero = { |m, rC, cC|
var sz = m.size, out = nil;
block { |br|
sz.do { |i|
if(rC[i].not) {
sz.do { |j|
if((m[i][j] == 0) && cC[j].not) {
out = [i, j]; br.value;
}
}
}
}
};
out
};
// Check a specific row for a starred zero.
~starRow = { |mask, r|
var sz = mask.size, c = -1;
sz.do { |j| if(mask[r][j] == 1) { c = j } };
c
};
// Find the smallest value that is not covered by any row or column
~minUnc = { |m, rC, cC|
var sz = m.size, mn = inf;
sz.do { |i|
if(rC[i].not) {
sz.do { |j| if(cC[j].not) { mn = min(mn, m[i][j]) } }
}
};
mn
};
// Adjust the matrix/array
~adjMat = { |m, rC, cC, mn|
var sz = m.size, r = m.deepCopy;
sz.do { |i|
sz.do { |j|
if(rC[i]) { r[i][j] = r[i][j] + mn };
if(cC[j].not) { r[i][j] = r[i][j] - mn };
}
};
r
};
// Iterates over the mask and for each row,
// if a starred zero is found, the column index is recorded.
// The resulting array maps each source (row) to a target (column)
~extract = { |mask|
var sz = mask.size, out = Array.fill(sz, -1);
sz.do { |i| sz.do { |j| if(mask[i][j] == 1) { out[i] = j } } };
out
};
// Main function
~hungarian = { |cost|
var sz = cost.size, m = ~colMin.(~rowMin.(cost));
var mk = ~starZero.(m), rC = Array.fill(sz, false), cC = Array.fill(sz, false);
var st = 3, dn = false, maxI = 100, it = 0;
while({ dn.not and: { it < maxI } }, {
it = it + 1;
switch(st,
3, {
rC = Array.fill(sz, false); cC = Array.fill(sz, false); st = 4;
},
4, {
var cCount;
cC = ~coverStar.(mk);
cCount = cC.count({ |item| item == true });
if(cCount >= sz) { dn = true } { st = 5 };
},
5, {
var z = ~findZero.(m, rC, cC);
if(z.isNil) {
st = 6;
} {
var row = z[0], col = z[1], sc = ~starRow.(mk, row);
mk[row][col] = 2;
if(sc >= 0) {
rC[row] = true; cC[sc] = false;
} {
st = 6;
};
};
},
6, {
var mn = ~minUnc.(m, rC, cC);
m = ~adjMat.(m, rC, cC, mn);
st = 4;
}
);
});
if(it >= maxI) {
"Maximum iterations without convergence.".warn;
};
~extract.(mk)
};
)
(
// Voice leading between two chords
~voiceLead = { |src, trg|
var cm = src.collect { |n1| trg.collect { |n2| (n1 - n2).abs } };
var asn = ~hungarian.(cm);
asn.collect { |t, s| [src[s], trg[t]] }
};
// Example
~chC = [60, 64, 67];
~chD = [62, 65, 69];
~cMat = ~chC.collect { |n1| ~chD.collect { |n2| (n1 - n2).abs } };
~cMat.do { |r| r.postln };
~vLead = ~voiceLead.(~chC, ~chD);
~vLead.do { |p|
var src = p[0], trg = p[1], int = trg - src;
("Note " ++ src ++ " → " ++ trg ++ " (interval: " ++ int ++ ")").postln;
};
)
2 Likes
Incredible! Thanks so very much, @smoge! It feels to me like it might be worth considering including this or something like it in the SC code. What do you think?
1 Like
Thanks. Glad it helps.
Well, it could be, or maybe a quark. Feel free to adapt it.
It would need some refactoring and organizing to turn into a class(es), methods, etc.
Also: tests, documentation, error handling, etc.
This implementation may fail with edge cases; maybe more work is needed to improve it.
Let me know how it performs with your materials
Hi @smoge, in fact I’ve just been testing it and I think there may be something amiss, unless I’ve overlooked something. The example on Wikipedia is this:
~costMatrix = [
[8, 4, 7],
[5, 2, 3],
[9, 4, 8]
];
With expected result [ 0, 1, -1] in the format your code provides, mapping source row to target column.
But your implementation outputs [ 1, 0, -1 ], which in any case seems wrong, since it would map all rows to the same column.
For comparison, here’s some brute-force code (so it solves in factorial time rather than in polynomial time) I wrote that gives the correct answer (in format [ [ 0, 0 ], [ 2, 1 ], [ 1, 2 ] ]):
+SequenceableCollection {
assignWorkersForLowestCost{
var minCost = inf;
var minConfiguration;
var instrumentIndices = Array.fill(this.size, {|i| i});
this.size.factorial.do{|i|
var cost = instrumentIndices.permute(i).collect({|item,j| this[item][j]}).sum;
if(cost<minCost,{
minCost = cost;
minConfiguration = instrumentIndices.permute(i);
});
};
^minConfiguration.collect({|item, i| [item, i] });
}
}
For your voice-leading example, CMaj triad —> DMaj triad, my brute-force code gives the same result your code does:
(
~costMatrix = [
[ 2, 5, 9 ],
[ 2, 1, 5 ],
[ 5, 2, 2 ]
];
~costMatrix.assignWorkersForLowestCost //(my brute force approach!)
)
with result: [ [ 0, 0 ], [ 1, 1 ], [ 2, 2 ] ]
Yes, there are some limitations.
One of them is here; this only subtracts when the minimum is positive, which is a potential issue:
if(mn > 0) { sz.do { |i| r[i][j] = r[i][j] - mn } };
I can try to review it sometime soon.
1 Like
Your “brute force” code can be pretty handy for debugging the “Hungarian” optimization.
EDIT: I believe the discrepancy can also arise because state 5 is incomplete. It does not attempt to find an augmenting path, so it may not always determine the optimal assignments for specific matrices. Maybe another hint on how to improve it. Although this enhancement is not trivial, it will require additional logic and many lines of code.
1 Like
Let me know if you have a moment to dig back in, @smoge, I’d be very interested! FYI, I had similar issues with a Python implementation that didn’t use libraries. On the other hand, the known Python libraries (like scipy’s linear_sum_assignment) do it perfectly, of course. So it’s clearly not trivial.
1 Like
Hi again @smoge, I just did a test of the brute-force method with larger cost matrices, and as expected, it’s brutal: about 6 seconds for a 10x10 matrix, which is a size that’s not unrealistic at all when you’re talking about assigning instruments to notes. If you’re inclined to take another look at your code, I’d be super grateful!
PS— in the mean time, I’m sending OSC to the scipy’s linear_sum_assignment library and back, and that’s working fine. The speed gains over brute force are pretty spectacular. Here’s code for that if anyone needs the Munkres in SC, like I do right now:
Python:
from pythonosc.dispatcher import Dispatcher
from pythonosc.osc_server import ThreadingOSCUDPServer
from pythonosc.udp_client import SimpleUDPClient
import numpy as np
from scipy.optimize import linear_sum_assignment
def hungarian_handler(address, *args):
size = int(args[0]) # explicitly ensure integer
flat_matrix = [float(x) for x in args[1:]] # explicitly ensure float list
cost_matrix = np.array(flat_matrix).reshape((size, size))
print(f"Received cost matrix:\n{cost_matrix}")
row_ind, col_ind = linear_sum_assignment(cost_matrix)
assignment = col_ind.tolist()
print(f"Optimal assignment (columns): {assignment}")
# Send results back to SuperCollider
client = SimpleUDPClient("127.0.0.1", 57120) # SuperCollider default port
client.send_message("/hungarian_result", assignment)
dispatcher = Dispatcher()
dispatcher.map("/hungarian_matrix", hungarian_handler)
ip = "127.0.0.1"
port = 5005
server = ThreadingOSCUDPServer((ip, port), dispatcher)
print(f"OSC server running on {ip}:{port}")
server.serve_forever()
and SC:
(
// Receive the results from Python
OSCdef(\hungarianResults, { |msg|
var results = msg[1..];
"Optimal Assignment from Python: %\n".format(results).postln;
}, '/hungarian_result');
//wikipedia example
~costMatrix = [
[8, 4, 7],
[5, 2, 3],
[9, 4, 8]
];
//expected: [0,2,1]
//send to Python:
~addr = NetAddr("127.0.0.1", 5005);
~addr.sendMsg("/hungarian_matrix", ~costMatrix.size, *~costMatrix.flat);
)
Finally, for anyone interested, here’s a very nice, well-commented implementation in Python:
1 Like
Hey, my man,
I borrowed the missing parts from Scipy (especially the path construction algorithm ). Can you test this version and compare it to the brute-force results?
(
~rowMin = { |m|
m.collect { |r|
var mn;
mn = r.minItem;
r.collect { |v| v - mn }
}
};
~colMin = { |m|
var sz, r;
sz = m[0].size;
r = m.deepCopy;
sz.do { |j|
var mn;
mn = inf;
m.size.do { |i| mn = min(mn, r[i][j]) };
if(mn > 0) { m.size.do { |i| r[i][j] = r[i][j] - mn } };
};
r
};
~starZero = { |m|
var n, cols, mask, rC, cC;
n = m.size;
cols = m[0].size;
mask = Array.fill(n, { Array.fill(cols, 0) });
rC = Array.fill(n, false);
cC = Array.fill(cols, false);
n.do { |i|
cols.do { |j|
if((m[i][j] == 0) && rC[i].not && cC[j].not) {
mask[i][j] = 1; rC[i] = true; cC[j] = true;
}
}
};
mask
};
~coverStar = { |mask|
var n, cols, cC;
n = mask.size;
cols = mask[0].size;
cC = Array.fill(cols, false);
cols.do { |j|
n.do { |i|
if(mask[i][j] == 1) { cC[j] = true }
}
};
cC
};
~findZero = { |m, rC, cC|
var n, cols, out;
n = m.size;
cols = m[0].size;
out = nil;
block { |br|
n.do { |i|
if(rC[i].not) {
cols.do { |j|
if((m[i][j] == 0) && cC[j].not) {
out = [i, j]; br.value;
}
}
}
}
};
out
};
~primeZero = { |mask, r, c|
var cols, starCol;
cols = mask[0].size;
starCol = -1;
mask[r][c] = 2;
block { |break| cols.do { |j|
if(mask[r][j] == 1) {
starCol = j;
break.value;
}
}};
starCol
};
~starRow = { |mask, r|
var cols, c;
cols = mask[0].size;
c = -1;
cols.do { |j| if(mask[r][j] == 1) { c = j } };
c
};
~starCol = { |mask, c|
var n, r;
n = mask.size;
r = -1;
n.do { |i| if(mask[i][c] == 1) { r = i } };
r
};
~minUnc = { |m, rC, cC|
var n, cols, mn;
n = m.size;
cols = m[0].size;
mn = inf;
n.do { |i|
if(rC[i].not) {
cols.do { |j|
if(cC[j].not) {
mn = min(mn, m[i][j])
}
}
}
};
mn
};
~constructPath = { |mask, r0, c0|
var n, cols, path, count, done, r, c;
n = mask.size;
cols = mask[0].size;
path = List.new;
count = 0;
done = false;
r = r0;
c = c0;
path.add([r, c]);
while { done.not } {
var starRow;
starRow = ~starCol.(mask, c);
if(starRow < 0) {
done = true;
} {
var primeCol;
r = starRow;
path.add([r, c]);
primeCol = -1;
block { |break| cols.do { |j|
if(mask[r][j] == 2) {
primeCol = j;
break.value;
}
}};
if(primeCol < 0) {
"No primed zero found in row".warn;
done = true;
} {
c = primeCol;
path.add([r, c]);
};
};
count = count + 1;
if(count > (n * cols)) {
"Maximum path length exceeded.".warn;
done = true;
};
};
path.asArray
};
~augmentPath = { |mask, path|
path.do { |pos|
var r, c;
r = pos[0];
c = pos[1];
if(mask[r][c] == 1) {
mask[r][c] = 0;
} {
mask[r][c] = 1;
};
};
mask.size.do { |i|
mask[0].size.do { |j|
if(mask[i][j] == 2) { mask[i][j] = 0 };
};
};
mask
};
~adjMat = { |m, rC, cC, mn|
var r;
r = m.deepCopy;
m.size.do { |i|
m[0].size.do { |j|
if(rC[i]) { r[i][j] = r[i][j] + mn };
if(cC[j].not) { r[i][j] = r[i][j] - mn };
}
};
r
};
~ensureShape = { |cost|
var m, transposed;
m = cost.deepCopy;
transposed = false;
if(m.size > m[0].size) {
var tCost;
tCost = Array.fill(m[0].size, { Array.fill(m.size, 0) });
m.size.do { |i|
m[0].size.do { |j|
tCost[j][i] = m[i][j];
}
};
m = tCost;
transposed = true;
};
[m, transposed]
};
~extract = { |mask, transposed|
var n, cols, rowInd, colInd;
n = mask.size;
cols = mask[0].size;
rowInd = List.new;
colInd = List.new;
n.do { |i|
cols.do { |j|
if(mask[i][j] == 1) {
if(transposed) {
rowInd.add(j);
colInd.add(i);
} {
rowInd.add(i);
colInd.add(j);
}
}
}
};
[rowInd.asArray, colInd.asArray]
};
~hungarian = { |cost|
var costMat, transposed, mask, rC, cC;
var n, cols, done, maxIter, iter;
#costMat, transposed = ~ensureShape.(cost);
n = costMat.size;
cols = costMat[0].size;
rC = Array.fill(n, false);
cC = Array.fill(cols, false);
costMat = ~colMin.(~rowMin.(costMat));
mask = ~starZero.(costMat);
done = false;
maxIter = 100;
iter = 0;
while { done.not && (iter < maxIter) } {
var starCount;
iter = iter + 1;
cC = ~coverStar.(mask);
starCount = cC.count({ |item| item == true });
if(starCount >= n) {
done = true;
} {
var z = ~findZero.(costMat, rC, cC);
if(z.isNil) {
var minVal;
minVal = ~minUnc.(costMat, rC, cC);
costMat = ~adjMat.(costMat, rC, cC, minVal);
} {
var row, col, starCol;
row = z[0];
col = z[1];
starCol = ~primeZero.(mask, row, col);
if(starCol < 0) {
var path;
path = ~constructPath.(mask, row, col);
mask = ~augmentPath.(mask, path);
rC = Array.fill(n, false);
cC = Array.fill(cols, false);
} {
rC[row] = true;
cC[starCol] = false;
};
};
};
};
if(iter >= maxIter) {
"Maximum iterations reached without convergence.".warn;
};
~extract.(mask, transposed)
};
~voiceLead = { |src, trg|
var costMatrix, rowInd, colInd, mapping;
costMatrix = src.collect { |n1| trg.collect { |n2| (n1 - n2).abs } };
mapping = List.new;
#rowInd, colInd = ~hungarian.(costMatrix);
rowInd.size.do { |i|
mapping.add([src[rowInd[i]], trg[colInd[i]]]);
};
mapping.asArray
};
)
(
~chordC = [60, 64, 67];
~chordF = [65, 69, 72];
~voiceLeading = ~voiceLead.(~chordC, ~chordF);
~voiceLeading.do { |pair, i|
var src, trg, interval;
src = pair[0];
trg = pair[1];
interval = trg - src;
("Voice " ++ (i+1) ++ ": " ++ src ++ " → " ++ trg ++ " (interval: " ++ interval ++ " semitones)").postln;
};
~chord1 = [60, 64, 67, 72];
~chord2 = [65, 69, 72];
~voiceLeading2 = ~voiceLead.(~chord1, ~chord2);
~voiceLeading2.do { |pair, i|
var src, trg, interval;
src = pair[0];
trg = pair[1];
interval = trg - src;
("Voice " ++ (i+1) ++ ": " ++ src ++ " → " ++ trg ++ " (interval: " ++ interval ++ " semitones)").postln;
};
)
Hey @smoge,
I’m not quite sure what I’m looking at in the results — the ~hungarian() function now seems to output two different mapping options, including in your voice-leading example. And when I run
~costMatrix = [
[8, 4, 7],
[5, 2, 3],
[9, 4, 8]
];
~hungarian.(~costMatrix)
with your new code, I get [ [ 0, 1, 2 ], [ 1, 0, 2 ] ] as a result. Neither one of the 2 options is the expected result of [0,1,-1].
1 Like
Yes, I changed some things:
costMatrix = [ [cost_1_1, cost_1_2, ...], [cost_2_1, cost_2_2, ...], ... ];
#rowIndices, colIndices = ~hungarian.(costMatrix);
The example gives me this:
Voice 1: 60 → 65 (interval: 5 semitones)
Voice 2: 67 → 69 (interval: 2 semitones)
Voice 1: 60 → 65 (interval: 5 semitones)
Voice 2: 67 → 69 (interval: 2 semitones)
Voice 3: 72 → 72 (interval: 0 semitones)
-> [[60, 65], [67, 69], [72, 72]]
The result [ [ 0, 1, 2 ], [ 1, 0, 2 ] ] means:
- Row 0 is assigned to column 1
- Row 1 is assigned to column 0
- Row 2 is assigned to column 2
EDIT: Later I will take a look if things are correct
Thanks much @smoge! I thought that might be what the new output format was doing. Unfortunately that means it’s not giving the right result (although it’s better than before when it was pointing all three workers to the same job). In the new notation, the expected result for the wikipedia example is: [ [ 0, 1, 2 ], [ 0, 2, 1 ] ].
1 Like
Ok, thanks
I will then try to find the bug. My suspects are ~constructPath, ~findZero, or something in the ordering operations in the main function.
1 Like
Wish I could help. Up to my ears for a big deadline. Could help investigate in April!
A couple years ago I did a port of Munkres to sclang. I also did a port of the “Auction” Algorithm, which also does linear assignment but finds a decent solution faster than Munkres (Munkres guarantees an optimal solution). I haven’t look at this in a while but I was using them at the time so should be good!
Let me know what you find…
4 Likes
Wow, this is fantastic. Tested it and it works beautifully. Thanks so very much! In my opinion, this implementation should be included in SC. A fundamental algorithm that is quite useful in music.
2 Likes
For what it’s worth - while SC is a proper, Turing Complete language, its strengths lie in domain specific audio processing: wiring UGens together for synthesis and signal processing, and creating sequences through Patterns and routines.
Decisions that need to be made at audio rate or control rate are best done in SC. Decisions/actions that take place at a slower rate (a “performance” rate), choosing instruments, pitches, changing effects, etc, can be farmed out to another system that communicates with SC via OSC. And if they are complicated decisions, I think they should be.
There are trade offs - and as always “there’s more than one way to do it”. Maybe your algorithm can be efficiently and precisely implemented in SC, but if it’s already implemented in a well debugged Python library, why not use that?
Similarly, you could implement graph theory algorithms or neural network algorithms in SC, but that seems like a huge waste of time. Again, it’s your time, and all power to you, but I think the answer to “why isn’t this in SC?” is, “that’s not what SC is for”.
2 Likes
Good point, well expressed 
. However, there are many things included in SC that don’t fit the paradigm you’re describing. So including a good port of a fundamental algorithm makes sense to me.
Ok, fair. But you used a bit of straw man rethoric there. I don’t think things happen in such complex and specialized case.
Another reason why, for example, Miller Puckett says that SC is “… a … programming language (but sort of…)” is because there is no culture of hacking and analysing algorithms. No one on the forum discusses how to make a quicksort more efficient in SC, for example. (I’m just using his example, but can be anything)