assignment.pi

import cp.
        
matrix_square_sub_total(M,M2) = M3 =>
    M3 = new_array(M.length, M.length),
    foreach(I in 1..M.length, J in 1..M.length)
        M3[I,J] = sum([pow(M[I,K]-M2[J,K],2) : K in 1..M[1].length])
    end.

% modeled after http://www.hakank.org/picat/assignment.pi
nearest(M1,M2) = Assignments =>
    Cost=matrix_square_sub_total(M1,M2),
    % get the dimension of the problem
    Rows = Cost.length,
    Cols = Cost[1].length,

    % decision variables: a 0..1 matrix
    X = new_array(Rows,Cols),
    X :: 0..1,

    % exacly one assignment per row, all rows must be assigned
    foreach(I in 1..Rows) sum([X[I,J] : J in 1..Cols]) #= 1 end,

    % zero or one assignments per column
    foreach(J in 1..Cols) sum([X[I,J] : I in 1..Rows]) #=< 1 end,

    % calculate TotalCost
    TotalCost #= sum([X[I,J]*Cost[I,J] : I in 1..Rows, J in 1..Cols]),
    
    %
    % search
    % 
    solve([$min(TotalCost)], X),

    %
    % get the assignments 
    %
    % note: this should be done _after_ labeling.
    %
    Assignments = [J : I in 1..Rows, J in 1..Cols, X[I,J] == 1],
    writeln(assignments=Assignments),
    writeln(total_cost=TotalCost),
    nl.

assignmentTest1 =>
    M1={{1,2,3},
        {100,200,300},
        {1000,2000,3000}},
        
    M2={
        {105,205,305},
        {4,6,8},
        {1005,2005,3010}},
    Assignments=nearest(M1,M2),
    %writeln(Assignments),
    %check results
    Assignments=[2,1,3].

assignmentTest2 =>
    M1={{1000},{50},{1},{100}},        
    M2={{1},{1000},{5},{50}},
    Assignments=nearest(M1,M2),
    % writeln(Assignments),
    % check results
    /*
        initially it seemed like the least cost may have been
           1000=1000,
           50=50,
           1=1,
           100=5
           
        Although that minimizes individual differences, the total
        squared diff is 9025
        (1000, 50, 1, 100) (+/@:*:@:-)  (1000,50,1,5)
        9025

        vs. 4525
           1000=1000,
           50=5,
           1=1,
           100=50

           (1000, 50, 1, 100) (+/@:*:@:-)  (1000,5,1,50)
           4525
        brute force in J:
        ({. /: sols) { space [ sols=: (1000, 50, 1, 100)&(+/@:*:@:-)"1 space [ space=:(i.!4) A. (1,5,50,1000)
    */
    Assignments=[2,3,1,4].


main=>
    assignmentTest1(),
    assignmentTest2(),
    true.