Preamble
I can't resist to answer, although the answer won't be short :-). Below is an attempt to systematically capitalize on properties which weren't used (or at least used systematically) in other answers.
I will use the same lists provided by the OP, for tests:
l1 = {5.7832, 30.4713, 74.887, 139.04, 222.932, 326.563}
l2 = {3.481, 9.2816, 15.7112, 27.8226, 45.0379, 67.222, 94.3526}
The algorithm
The main idea is that for many common binary operators, the order of the list resulting from applying the operator to a list and a single element of another list can be deduced from the order of the original list and the element in question, often without explicit sorting of the resulting list.
The algorithm itself will be based on lazy merging of lazy lists, rsulting from application of the binary operator in question to a list (say l1
) and each of the elements of the other list l2
. What we will do is to create a binary search tree,which holds the information nodes of the following form:
{{currentNumber, index1,index2},iterator}
where an iterator is a construct which supplies the next element of a resulting list on demand. What is important is that an iterator encodes the information about the ordering of the new sub-list.
What we will do is to merge the resulting sub-lists lazily, very analogously to what merge
function does for 2 lists in the mergesort algorithm, but now we have Length[l2]
sub-lists (represented by their iterators). Here is the outline of the general algorithm to get a single element with indices:
Create Length[l2]
iterators, each encoding a lazy version of the result of an application of our binary operator to a list l1
and i
-th element of l2
.
Construct an initial binary search tree by requesting each iterator to provide the corresponding first element and constructing nodes of the form outlined above. The ordering used for the binary seacrh tree construction can be also specified.
If we already picked enough elements, exit.
Pick the node with the smallest (according to our ordering) element, and save it in a temporary variable.
Delete this node from a tree.
Construct the new node by requesting the next element from the iterator of the smallest (deleted) node, if that iterator still has elements, and add this node to the tree.
Return the first part of the deleted node, {currentNumber, index1,index2}
Go to 3
Since insertion and deletion time for a (balanced) binary search tree is O(Log(m))
, where m
is the number of nodes (which is the current number of iterators, can be at most Length[l2]
), the complexity of the search alone will be of the order of n Log[m]
, where n
is the number of elements requested.
This can be very fast for small enough n
and large enough lists, but the main question is whether or not we can build the iterators efficiently. In general, this would involve copying and sorting all sub-lists, and therefore should have a complexity Length[l1]Log[Length[l1]]Length[l2]
. However, for some common cases (such as Times and
Plus`), the resulting ordering can be deduced from the initial ordering and the element of the second list in just constant time, and then we win big both speed-wise and memory-wise.
Mathematica impelentation
Binary search tree
For simplicity, I will use the unbalanced binary search tree implementation, which is a slightly generalized (to allow a custom ordering function) version of an excellent implementation due to Roman Maeder ("Computer Science With Mathematica", chapter 6). Here it is:
These are constructors and selectors for the tree nodes:
Clear[node, tree, leftTree, rightTree, emptyTree];
node[info_, left_, right_] := tree[info, left, right]
leftTree[tree[_, left_, _]] := left
rightTree[tree[_, _, right_]] := right
information[tree[info_, _, _]] := info
This is the insertion function:
Clear[insertTree];
insertTree[emptyTree, info_, Order_: Order, Key_: Identity] :=
node[info, emptyTree, emptyTree]
insertTree[tree_, info_, Order_: Order, Key_: Identity] /;
Order[Key[info], Key[information[tree]]] == 0 := tree;
insertTree[tree_, info_, Order_: Order, Key_: Identity] /;
Order[Key[info], Key[information[tree]]] > 0 :=
node[
information[tree],
insertTree[leftTree[tree], info, Order, Key],
rightTree[tree]
];
insertTree[tree_, info_, Order_: Order, Key_: Identity] :=
node[
information[tree],
leftTree[tree],
insertTree[rightTree[tree], info, Order, Key]
]
Node insertion is a relatively simple matter since one just needs to find the place where to insert it, and the way the tree is constructed automatically preserves the right order.
We won't need the node search functionality, so I skip it. This is the function to find the smallest node by straightforward tree descent:
Clear[smallestNode];
smallestNode[emptyTree] = Null
smallestNode[tree_] /; leftTree[tree] === emptyTree :=
information[tree]
smallestNode[tree_] := smallestNode[leftTree[tree]]
Here is the most conceptually difficult function to delete the node:
Clear[deleteTree];
deleteTree[tree_, key_, Order_: Order, Key_: Identity] /;
Order[Key[information[tree]], key] < 0 :=
node[
information[tree],
deleteTree[leftTree[tree], key, Order, Key],
rightTree[tree]
]
deleteTree[tree_, key_, Order_: Order, Key_: Identity] /;
Order[Key[information[tree]], key] > 0 :=
node[
information[tree],
leftTree[tree],
deleteTree[rightTree[tree], key, Order, Key]
]
deleteTree[tree_, key_, Order_: Order, Key_: Identity] /;
leftTree[tree] == emptyTree := rightTree[tree]
deleteTree[tree_, key_, Order_: Order, Key_: Identity] /;
rightTree[tree] == emptyTree := leftTree[tree]
deleteTree[tree_, key_, Order_: Order, Key_: Identity] :=
With[{nextinfo = smallestNode[rightTree[tree]]},
node[
nextinfo,
leftTree[tree],
deleteTree[rightTree[tree], Key[nextinfo], Order, Key]
]
]
Finally, here is the constructor to construct tree from a list of information items:
Clear[Tree];
Tree[l_List, Order_: Order, Key_: Identity] :=
Fold[Function[{t, i}, insertTree[t, i, Order, Key]], emptyTree, l]
Some examples:
(t = Tree[{5, 4, 7, 3, 2, 8, 6}]) // TreeForm
deleteTree[t, 5] // TreeForm
I may include some visualization function to plot a tree better than TreeForm
, in the future.
Iterators
Here is a generic iterator "data type":
Clear[getNext, iterator];
getNext[iterator[fun_]] := fun[];
iterator /: MakeBoxes[it_iterator, fmt_] :=
InterpretationBox[RowBox[{"«", "iterator", "»"}], it];
It basically tells that objects with heads iterator
contain a function embedded in them, and the end user should only know that applying to the getNext
, s/he gets the next element.
We now become specific and construct the iterator generator for iterators for Plus
binary operator:
Clear[makePlusInterator];
makePlusInterator[l1_, elem_] :=
With[{len = Length[l1]},
Module[{i = 0},
iterator[Function[If[i < len, l1[[++i]] + elem, None]]]
]
];
makePlusInterator[l1_, l2_List] :=
Map[makePlusInterator[l1, #] &, l2];
Our iterators so constructed are closures, embedding addition as a part of producing an element on demand. Here we encode our knowledge of the resulting ordering (which, in the case of Plus
, is the same as the original one).
Here is an example:
iter = makePlusInterator[l2,First@l1]
(* « iterator » *)
Table[getNext[iter],{6}]
(* {9.2642,15.0648,21.4944,33.6058,50.8211,73.0052} *)
Note that iterators themselves do not carry the index information (perhaps it would be more logical if they would).
Bringing it together
To combine the components I just described and actually implement the algorithm, we need two functions. First is a helper function which implements a part of the step 6:
Clear[getNextInfo];
getNextInfo[{{next_, fstind_, secind_}, iter_}] :=
{{getNext[iter], fstind, secind + 1}, iter};
It takes the information stored in a node and prepares the information for the new node to be inserted into the tree. Here is the main function:
Clear[makeResultinTreeInterator];
makeResultinTreeInterator[iterators : {__iterator}, Order_: Order] :=
Module[{
tr =
Tree[
MapIndexed[{{getNext[#], First@#2, 1}, #} &, iterators],
Order,
First
]
},
iterator[
Function[
Module[{smallest = smallestNode[tr], next},
tr = deleteTree[tr, smallest];
next = getNextInfo[smallest];
If[First@next =!= None,
tr = insertTree[tr, next]
];
First@smallest]
]
]
];
Note that the function itself returns an iterator, so our resulting list is lazy as well.
Examples
I will now illustrate this:
resiter = makeResultinTreeInterator[makePlusInterator[l2, l1], Order]
Table[getNext[resiter],{6}]
Table[getNext[resiter],{6}]
Table[getNext[resiter],{6}]
(*
{{9.2642,1,1},{15.0648,1,2},{21.4944,1,3},{33.6058,1,4},{33.9523,2,1},{39.7529,2,2}}
{{46.1825,2,3},{50.8211,1,5},{58.2939,2,4},{73.0052,1,6},{75.5092,2,5},{78.368,3,1}}
{{84.1686,3,2},{90.5982,3,3},{97.6933,2,6},{100.136,1,7},{102.71,3,4},{119.925,3,5}}
*)
So far so good. Now, we pick some larger lists:
large1 = Sort@RandomReal[{-100, 100}, 1000];
large2 = Sort@RandomReal[{-100, 100}, 100];
randomizedOrderLarge2 = RandomSample[large2];
resiterLarge =
makeResultinTreeInterator[
makePlusInterator[large1 , randomizedOrderLarge2], Order];
Here I had to randomize the second sample because the implementation of the binary search tree I use is for unbalanced trees, which would degenerate to a line for an ordered second list. This somewhat violates the initial condition that both lists are ordered, but I only use this for benchmarks, and the resulting list won't change in any case.
We now test by fetching some larger number of elements:
Table[getNext[resiterLarge],{500}]//Short//AbsoluteTiming
Table[getNext[resiterLarge],{1000}]//Short//AbsoluteTiming
Table[getNext[resiterLarge],{2000}]//Short//AbsoluteTiming
(*
{0.3925781,{{-191.248,33,1},<<498>>,{-172.676,27,43}}}
{0.8896485,{{-172.665,6,53},<<998>>,{-161.276,29,162}}}
{1.9609375,{{-161.251,56,14},<<1998>>,{-144.795,20,134}}}
*)
We see that the timing linearly depends on the number of elements, but our code is really slow.
Java implementation
Now that we have our slow Mathematica prototype, we can port it to Java in the hopes that this will run (much) faster. To work with this part, you will need the Java reloader (which I modified to recognize interfaces, so if you have the old one you should change to use the recent version).
While what follows is conceptually a direct port of the Mathematica part, there will be significant differences in details. In particular, I will use the Java TreeSet
library class, so that I don't need to implement the binary search tree from scratch.
So, run the Java reloader code first.
General interface for customization
The following interface should be implemented for each binary operator and ordering, to customize the algorithm (execute this code):
JCompileLoad@
"package orderedPairsFinal;
import java.util.Comparator;
import java.util.Iterator;
public interface Customizer {
public Iterator<Double> getIterator(double[] list, double elem);
public Comparator<Double> getComparator();
}"
It basically declares that we have a way to compare doubles and construct iterators (I stick to doubles here).
Customization for Plus
This implements the above interface for Plus
:
JCompileLoad@"
package orderedPairsFinal;
import java.util.Comparator;
import java.util.Iterator;
public class PlusCustomizer implements Customizer {
@Override
public Comparator<Double> getComparator() {
return new Comparator<Double>() {
public int compare(Double arg0, Double arg1) {
return Double.compare(arg0, arg1);
}
};
}
@Override
public Iterator<Double> getIterator(final double[] list,
final double elem) {
return new Iterator<Double>() {
private int ctr = 0;
public boolean hasNext() {
return ctr < list.length - 1;
}
public Double next() {
return list[ctr++] + elem;
}
public void remove() { }
};
}
}"
Tree node class
The following class implements the tree node abstraction.
JCompileLoad@"
package orderedPairsFinal;
import java.util.Comparator;
import java.util.Iterator;
public class NodeInfo implements Iterator<NodeInfo>,Comparable<NodeInfo> {
private Iterator<Double> iter;
private Comparator<NodeInfo> comp;
public Double current;
public final int fstIndex;
public int secIndex;
public NodeInfo(Iterator<Double> iter,Comparator<NodeInfo> comp,
int index, int secIndex) {
this.fstIndex = index;
this.secIndex = secIndex;
this.iter = iter;
this.comp = comp;
current = iter.next();//Note: no checks performed here
}
@Override
public NodeInfo next() {
return new NodeInfo(iter,comp,fstIndex,secIndex+1);
}
@Override
public boolean hasNext() {
return iter.hasNext();
}
public void remove() { }
@Override
public int compareTo(NodeInfo arg0) {
return comp.compare(this, arg0);
}
}"
The main class
This is the class which implements our algorithm.
JCompileLoad@"package orderedPairsFinal;
import java.util.Comparator;
import java.util.TreeSet;
public class OrderedPairs {
private double[] first;
private double[] second;
private Comparator<NodeInfo> comp;
private TreeSet<NodeInfo> tree = new TreeSet<NodeInfo>(comp);
public OrderedPairs(Customizer cust, double[] list1,
double[] list2) {
first = list1.length>list2.length?list1:list2;
second = list1.length>list2.length?list2:list1;
final Comparator<Double> dcomp = cust.getComparator();
comp =
new Comparator<NodeInfo>() {
public int compare(NodeInfo o1, NodeInfo o2) {
return dcomp.compare(o1.current, o2.current);
}
};
for(int i=0;i<second.length;i++){
tree.add(new NodeInfo(cust.getIterator(first, second[i]),comp,i,0));
}
}
public double[] getNext(){
if(tree.isEmpty()){
return null;
}
NodeInfo fst = tree.first();
tree.remove(fst);
if(fst.hasNext()){
tree.add(fst.next());
}
return new double[]{fst.current, fst.fstIndex+1,fst.secIndex+1};
}
public double[][] getNext(int num){
double[][] aux = new double[num][3];
int ctr = 0;
while(ctr < num){
double[] next = getNext();
if(next == null) break;
aux[ctr++] = next;
}
double[][] res = new double[ctr][3];
System.arraycopy(aux, 0, res, 0, ctr);
return res;
}
} "
The methods of this class and the way it works require some explanation. The class takes a customizer and two arrays of doubles in constructor, and creates iterators and the tree also in constructor. It has the method getNext
, which is overloaded to have 2 forms. You can either request a single element (with indices), or do that in a bulk. The latter possibility is particularly relevant for uses from Mathematica, to avoid the killing overhead of many individual method calls in JLink.
Examples
To use this, call:
cust = JavaNew["orderedPairsFinal.PlusCustomizer"];
op = JavaNew["orderedPairsFinal.OrderedPairs", cust, l1, l2];
We can now use it:
op@getNext[6]
op@getNext[6]
(*
{{9.2642,1.,1.},{15.0648,1.,2.},{21.4944,1.,3.},{33.6058,1.,4.},
{33.9523,2.,1.},{39.7529,2.,2.}}
{{46.1825,2.,3.},{50.8211,1.,5.},{58.2939,2.,4.},{73.0052,1.,6.},
{75.5092,2.,5.},{78.368,3.,1.}}
*)
Now, the larger benchmarks:
large1 = Sort@RandomReal[{-100, 100}, 1000];
large2 = Sort@RandomReal[{-100, 100}, 1000];
opL = JavaNew["orderedPairsFinal.OrderedPairs", cust, large1, large2];
and we can test:
opL@getNext[50000]//Short//AbsoluteTiming
opL@getNext[50000]//Short//AbsoluteTiming
(*
{0.0732421,{{-199.219,1.,1.},<<49998>>,{-137.102,73.,243.}}}
{0.0419922,{{-137.102,17.,304.},<<49998>>,{-110.502,335.,108.}}}
*)
which is quite fast for lazy lists to my mind, and this time includes the data transfer!(possibly being dominated by it).
Conclusions
I presented the abstract algorithm which promises good overall performance and very good performance for practically important special cases. I have also presented two implementations of it, using Mathematica and Java.
While Mathematica implementation has the right algorithmic behavior, it is very slow. This is one of the cases where the symbolic overhead of full Mathematica evaluator really gets in the way. The data structures involved and the structuring of code which allows to separate customizable parts from general parts (iterators, lazy lists, binary search tree) here did not allow to easily Compile
any part of the algorithm, so there are little chances to improve its speed within Mathematica.
Implementation in Java allowed to use powerful built-in data structures and idioms which in Java do not cost such an overhead as in Mathematica. As a result, the Java version is thousands times faster.
In both cases, I could also include the customizations for other common binary operators such as Times
, but did not do so simply because this post is already excessively long.