Lapack

General
1

Suppose I wanted to add a lapack binding for swipl. I don’t think such a thing exists, or does it? What would be a good strategy? I guess I would put it on top of either clapack or lapack++ (the c++ binding) to avoid having to interface with fortran. I am not sure, but I guess lapack++ would fit my needs even if some functions are missing. The matrix objects would be swipl “blobs” (?), with some accessor functions for row or column vectors etc.

Of course, there already exists a binding to R, which I could extend for my needs.

2

Roughly, I think there are three routes:

  • As you suggest, use it through R
  • Use it through the new Janus Python interface
  • Add a C(++) interface.

The latter is surely the most elaborate. Indeed, you represent external objects as blobs, either by packing a pointer or the data itself. In the end this route will give the best performance. The Janus approach could be the most attractive. Probably almost no work and pretty ok performance.

3

The annoying thing with R is that the type of the return values is so unpredictable. I have no experience with Janus, though, and I don’t really know python.

I guess I’ll try a little proof-of-concept with C++ and see how much work it is in terms of code lines.

4

You might be surprised how easy it is if there is a good Python binding. I’m not sure about that.

5

If you just need matrix ops, numpy is more or less the same as R matrix ops. All the usual slicing and indexing is supported. You can use numpy through Janus.

6

I am almost certain LAPACK and BLAS are only transitive dependencies of Numpy or Scipy or similar Python libraries

8

I’m sure you can :slight_smile: Do you have to have nice examples to share, for example in your blog?

9

I’ve nothing handy but I’ll rustle something up for you in the next few days :smiling_face:

1 Like
10

IMHO, the C++ interface (SWI-cpp2.h) is a lot easier than the C interface (the C interface requires correct code to do a lot of error checking; this is done automatically by the C++ interface). There’s also a new PlBlob class that simplifies implementing blobs.

But there’s also a bit of a learning curve, and the documentation could use some improvement. :frowning:

There’s also scipy, which seems to have an interface to LAPACK.

11

I tried the PlBlob sample code,

#include <SWI-cpp2.h>

PREDICATE(add, 3)
{ return A3.unify_integer(A1.as_long() + A2.as_long());
}


struct MyBlob;

static PL_blob_t my_blob = PL_BLOB_DEFINITION(MyBlob, "my_blob");

struct MyBlob : public PlBlob
{ std::unique_ptr<MyConnection> connection;
  std::string name_; // Used for error terms

etc., but I get an error when compiling it (ubuntu):

matthias@rechner:~/plblas$ ~/swipl-devel/build/src/swipl-ld -c plblas.cpp 
plblas.cpp:10:53: error: expected primary-expression before ‘,’ token
   10 | static PL_blob_t my_blob = PL_BLOB_DEFINITION(MyBlob, "my_blob");
      |                                                     ^
plblas.cpp:10:28: error: ‘PL_BLOB_DEFINITION’ was not declared in this scope; did you mean ‘PL_BLOB_VERSION’?
   10 | static PL_blob_t my_blob = PL_BLOB_DEFINITION(MyBlob, "my_blob");
      |                            ^~~~~~~~~~~~~~~~~~
      |                            PL_BLOB_VERSION
plblas.cpp:13:1: error: expected class-name before ‘{’ token
   13 | { std::unique_ptr<MyConnection> connection;
      | ^

Am I missing some include file?

12

I did a copy&paste from test_cpp.cpp and it compiled fine. Perhaps you’re missing the definition of MyConnection?

(It appears that the documentation is slightly behind the test code; I’ll update the documentation)

Here’s what worked for me, using the command swipl-ld:

sample blob code 
#include <SWI-cpp2.h>

PREDICATE(add, 3)
{ return A3.unify_integer(A1.as_long() + A2.as_long());
}

struct MyConnection
{ std::string name;

  explicit MyConnection() { }
  explicit MyConnection(const std::string& _name)
    : name(_name) { }

  ~MyConnection() { }

  bool open()
  { if ( name == "FAIL" ) // for testing error handling
      return false;
    return true;
  }

  bool close() noexcept
  { if ( name == "FAIL_close" ) // for testing error handling
      return false;
    return true;
  }
};


// The following code is taken from
// pl2cpp2.doc \subsubsection{Sample PlBlob code}
// with some minor changes for testing

struct MyBlob;

static PL_blob_t my_blob = PL_BLOB_DEFINITION(MyBlob, "my_blob");

struct MyBlob : public PlBlob
{ std::unique_ptr<MyConnection> connection;
  std::string name_; // Used for error terms

  explicit MyBlob()
    : PlBlob(&my_blob) { }

  explicit MyBlob(const std::string& connection_name)
    : PlBlob(&my_blob),
      connection(std::make_unique<MyConnection>(connection_name)),
      name_(connection_name)
  { if ( !connection->open() )
      throw MyBlobError("my_blob_open_error");
  }

  PL_BLOB_SIZE

  ~MyBlob() noexcept
  { if ( !close() )
      Sdprintf("***ERROR: Close MyBlob failed: %s\n", name_.c_str()); // Can't use PL_warning()
  }

  bool close() noexcept
  { if ( !connection )
      return true;
    bool rc = connection->close();
    connection.reset(); // Can be omitted, leaving deletion to ~MyBlob()
    return rc;
  }

  PlException MyBlobError(const char* error) const
  { return PlGeneralError(PlCompound(error, PlTermv(symbol_term())));
  }

  int compare_fields(const PlBlob* _b_data) const override
  { auto b_data = static_cast<const MyBlob*>(_b_data); // See note about cast
    return name_.compare(b_data->name_);
  }

  bool write_fields(IOSTREAM *s, int flags) const override
  { PlStream strm(s);

    strm.printf(",name=%s", name_.c_str());
    if ( !connection )
      strm.printf(",closed");
    return true;
  }
};

// %! create_my_blob(+Name: atom, -MyBlob) is semidet.
PREDICATE(create_my_blob, 2)
{ // Allocating the blob uses std::unique_ptr<MyBlob> so that it'll be
  // deleted if an error happens - the auto-deletion is disabled by
  // ref.release() before returning success.

  auto ref = std::unique_ptr<PlBlob>(new MyBlob(A1.as_atom().as_string()));
  return A2.unify_blob(&ref);
}

// %! close_my_blob(+MyBlob) is det.
// % Close the connection, silently succeeding if is already
// % closed; throw an exception if something goes wrong.
PREDICATE(close_my_blob, 1)
{ auto ref = PlBlobV<MyBlob>::cast_ex(A1, my_blob);
  if ( !ref->close() )
    throw ref->MyBlobError("my_blob_close_error");
  return true;
}

int main(void) { return 0; }
13

Problem solved. I think PL_BLOB_DEFINITION appeared lately, it is not yet in 9.0.4.

For the records: I compile the current swipl-devel (in subfolder build, say). I don’t install it yet. A working swipl is found in swipl-devel/build/src/swipl, also swipl-ld. The latter does not seem to adjust the include path properly. Edit: it uses the /usr/bin/swipl found in the path. Probably a feature, not necessarily a bug.

matthias@rechner:~/plblas$ ~/swipl-devel/build/src/swipl-ld -c plblas.cpp -v
	eval `swipl --dump-runtime-variables`
		PLBASE="/usr/lib/swi-prolog" <<<<<<<<<
		PLARCH="x86_64-linux"
		PLSOEXT="so"
		PLLIBDIR="/usr/lib/swi-prolog/lib/x86_64-linux"
		PLLIB="-lswipl"
		PLTHREADS="yes"
	/usr/bin/c++ -c -D_REENTRANT -D__SWI_PROLOG__ -D__SWI_EMBEDDED__ ***-I/usr/lib/swi-prolog/include*** -o plblas.o plblas.cpp
14

You need a fairly recent “development” version. What do you get from this?
swipl --version

15

Assuming you have a python3 environment, get numpy by typing pip3 install numpy at the command line, or python3 -m pip install numpy or perhaps py -m pip install numpy if under Windows. Here is some code:

% Doesn't try to convert complex Python objects to Prolog.
raw(Cmd,Out) :- py_call(Cmd,Out,[py_object(true)]).

% Tries to convert Python objects to Prolog.
cmd(Cmd,Out) :- py_call(Cmd,Out).

% Multiply two vectors.
example1 :-
    X = [1,2,3,4],
    Y = [2,3,4,5],
    raw(numpy:array(X),V1),
    raw(numpy:array(Y),V2),
    raw(numpy:multiply(V1,V2),V3),
    cmd(numpy:ndarray:tolist(V3), Out),   
    writeln(Out).

% Square a matrix.
example2 :-
    X = [[1,2,3],
	 [3,5,3],
	 [0,2,1]],
    raw(numpy:array(X),M),
    raw(numpy:matmul(X,X),Square),
    cmd(numpy:ndarray:tolist(Square),Out),
    writeln(Out).

It works as advertised by I’d note that getting something wrong (e.g. a typo) in py_call results in the REPL going into an infinite loop, but I’m sure that will be ironed out in due course.

There is an API reference for the vast array of function calls available to you here:

https://numpy.org/doc/stable/reference/index.html

If you want to put together a specific example I may be able to help!

NB: The NumPy linear algebra functions rely on BLAS and LAPACK to provide efficient low level implementations of standard linear algebra algorithms. Those libraries may be provided by NumPy itself using C versions of a subset of their reference implementations but, when possible, highly optimized libraries that take advantage of specialized processor functionality are preferred. Examples of such libraries are OpenBLAS, MKL ™, and ATLAS. Because those libraries are multithreaded and processor dependent, environmental variables and external packages such as threadpoolctl may be needed to control the number of threads or specify the processor architecture.

1 Like
16

Impressive, indeed. I added a use_module-directive, and removed a line in Example 2 that was copied from Example 1 (?).

use_module(library(janus)).

% Does not try to convert complex Python objects to Prolog.
raw(Cmd,Out) :-
    py_call(Cmd,Out,[py_object(true)]).
% Tries to convert Python objects to Prolog.
cmd(Cmd,Out) :- py_call(Cmd,Out).
% Square a matrix.
example2 :-
    X = [[1,2,3], [3,5,3], [0,2,1]],
    raw(numpy:transpose(X), XT),
    raw(numpy:matmul(XT, X),Square),
    cmd(numpy:ndarray:tolist(Square),Out),
    writeln(Out).

Can I get rid of the namespace numpy somehow? And is there a way to have two functions in the same expression? E.g. matmul(t(X), X)

1 Like
17

I had seen this syntax elsewhere using the janus interface, I think I can understand better what it means. The fact that “it tries” to convert Python objects means that the conversion may not succeed? And what one does get as an output in such cases? An error message? Thanks

PS: I mean: the call that doesn’t try to convert should always succeed, right?

18

Thanks!

Sure, I tried a few things though, but all works nicely. Any example?

Well, rather than raw/2 and cmd/2, you could define numpy/2, etc. You can pass the results of a Python function call into an argument using eval(Expr), e.g.

ex1 :-
    X = [1,2,3,4],
    Y = [2,3,4,5],
    py_call(numpy:ndarray:tolist(eval(numpy:multiply(eval(numpy:array(X)),
                                                     eval(numpy:array(Y))))),
            Out),
    writeln(Out).

If you want to use this for serious Prolog based linear algebra you probably want to define a DSL that concisely expresses your intend and interpret or compile this to janus calls. Maybe you can implement APL :slight_smile:

Conversion from Python to Prolog always succeeds (except for resource errors or Python exceptions in case an iterator is translated into a list). Without py_object(true), a number of Python objects are translated to natural counterparts in Prolog. Others produce, as with this option, a Python object reference.

19

I’m running version 9.1.20-9-g4344213cf under emacs. Starts like this:

?- py_call(zzz,_).
% Interative session; added `.` to Python `sys.path`
ERROR: Python 'ModuleNotFoundError':
ERROR:   No module named 'zzz'
ERROR: In:
ERROR:   [12] janus:py_call(zzz,_156)
ERROR:   [11] toplevel_call('<garbage_collected>') at /home/me/lib/swipl/boot/toplevel.pl:1317
^  Call: (6) forall(prolog:repl_loop_hook(end, _28), true) ?

I get an option regarding what to do, I pick “a” for abort. It then goes into an infinite loop, each segment of which looks like this:

ERROR: Arguments are not sufficiently instantiated
ERROR: In:
ERROR:    [9] _1322>=0
ERROR:    [8] prolog:repl_loop_hook(begin,_1356) at /home/me/lib/swipl/library/janus.pl:869
ERROR:    [7] forall(prolog:repl_loop_hook(begin,_1394),'$toplevel':true) at /home/me/lib/swipl/boot/apply.pl:52
ERROR:    [6] <meta call>
ERROR:    [5] sig_atomic('$toplevel':forall(...,true)) <foreign>
ERROR:    [4] setup_call_cleanup('$toplevel':forall(...,true),'$toplevel':'$query_loop'(0),'$toplevel':forall(...,true)) at /home/me/lib/swipl/boot/init.pl:679
20

Here is an example:

?- py_call(numpy:array([1,2,3]),X).
% Interative session; added `.` to Python `sys.path`
X = [<py_int64>(0x7f6963261950), <py_int64>(0x7f6963261d90), <py_int64>(0x7f69632618f0)].

To me this seems like a partial conversion, since the values inside of the list are not Prolog types. Meanwhile, X can’t be put back into any numpy function which expect an ndarray structures because its been converted to a list, so its practical to receive references instead and then convert at the end.

21

That has already been fixed shortly after it was introduced. Just bad luck :slight_smile: Upgrade and this should be fixed.

I saw this. It is a rather unfortunate result of the default translation. The numbers are instances of numpy.int64, which subclasses eventually from numpy.number, but this class hierarchy seems unrelated to Pythons own int class. As a result, Janus does not translate to integers, which would be nice. As is, the result is rather useless: it is neither useful for Python nor for the Prolog programmer.

Should we allow adding a method, e.g., __janus__() to any class that would then be used to obtain a value that is subsequently translated into a Prolog term? And, possibly, if __janus__() returns None, return an object reference? That would allow us extend libraries to behave the way we want.

Hmmm. It seems numpy.int64 is an immutable type and we cannot dynamically add methods to this. Some advice?