| 20 | |
| 21 | |
| 22 | int EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G, |
| 23 | double* alpha, double* beta, double *cost, uint64_t maxIter) { |
| 24 | // beware M and C are stored in row major C style!!! |
| 25 | |
| 26 | using namespace lemon; |
| 27 | uint64_t n, m, cur; |
| 28 | |
| 29 | typedef FullBipartiteDigraph Digraph; |
| 30 | DIGRAPH_TYPEDEFS(Digraph); |
| 31 | |
| 32 | // Get the number of non zero coordinates for r and c |
| 33 | n=0; |
| 34 | for (int i=0; i<n1; i++) { |
| 35 | double val=*(X+i); |
| 36 | if (val>0) { |
| 37 | n++; |
| 38 | }else if(val<0){ |
| 39 | return INFEASIBLE; |
| 40 | } |
| 41 | } |
| 42 | m=0; |
| 43 | for (int i=0; i<n2; i++) { |
| 44 | double val=*(Y+i); |
| 45 | if (val>0) { |
| 46 | m++; |
| 47 | }else if(val<0){ |
| 48 | return INFEASIBLE; |
| 49 | } |
| 50 | } |
| 51 | |
| 52 | // Define the graph |
| 53 | |
| 54 | std::vector<uint64_t> indI(n), indJ(m); |
| 55 | std::vector<double> weights1(n), weights2(m); |
| 56 | Digraph di(n, m); |
| 57 | NetworkSimplexSimple<Digraph,double,double, node_id_type> net(di, true, (int) (n + m), n * m, maxIter); |
| 58 | |
| 59 | // Set supply and demand, don't account for 0 values (faster) |
| 60 | |
| 61 | cur=0; |
| 62 | for (uint64_t i=0; i<n1; i++) { |
| 63 | double val=*(X+i); |
| 64 | if (val>0) { |
| 65 | weights1[ cur ] = val; |
| 66 | indI[cur++]=i; |
| 67 | } |
| 68 | } |
| 69 | |
| 70 | // Demand is actually negative supply... |
| 71 | |
| 72 | cur=0; |
| 73 | for (uint64_t i=0; i<n2; i++) { |
| 74 | double val=*(Y+i); |
| 75 | if (val>0) { |
| 76 | weights2[ cur ] = -val; |
| 77 | indJ[cur++]=i; |
| 78 | } |
| 79 | } |