Lapack

emiruz · December 17, 2023, 4:46pm

The basic object for numpy operations is the ndarray which has the tolist() function. The function converts the array (whatever its shape) into a list representation of basic Python types, which in turn converts just fine in Janus. So something like this works practically:

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), % <-- Convert as a last step.
    writeln(Out).

In this case at least, I don’t think that it is a problem, it just makes the py_object(true) more useful behaviour.

I personally think that py_object(true) is a better default. Say X has type T, when you convert X from Python to Prolog to Python, it may not have type T, and it may not be possible to pass it to other functions. That sort of thing is likely to happen with libraries like SpaCy, pytorch, tensorflow, and so on, which use complex types. I guess its going to happen whenever the response is a class … It happens in numpy too since an ndarray is iterable and will end up converted to a list… but then if you want to pass it to some function which expects to be able to run X.shape on its argument (ndarray has it, lists do not), the code will fail.

jan · December 17, 2023, 5:24pm

I had to fight hard to get py_object(true) into the specs for Janus The original idea is to write Python functions such that all the plumbing is done in Python and you merely call the function you defined with Prolog term on which you wanted some Python code to work, returning a Python object in Prolog-friendly format, so Prolog rarely sees a Python object reference. There is a lot to say for that, in particular in XSB as XSB does not support the blob interface which allows SWI-Prolog to represent handles of external objects safely and garbage collect them. So, a Python object in XSB needs to call py_free/1 on every object reference it sees exactly once. Failure to do so results in a memory leak. In SWI-Prolog you can call py_free/1, but you can also leave the object to GC. Prolog atom-gc will tell Python that the object has become unreachable in Prolog, so Python can do its GC on it.

That allows for safe management of Python objects from Prolog. Indeed, in this scenario returning object references by default, possibly except for base classes such as int makes sense. How much you want to do in Prolog and how much in Python is an open question. Note that the SWi-Prolog version also provides py_module/2 which allows you to keep your Python code close to the Prolog code. Note that, when using Prolog from C/C++/Python/JavaScript/Rust/… my advice is also to define a Prolog predicate that does all the plumbing and exchanges data in a format that is easy to handle on the other side.

I advocated for a bit more power in calling Python from Prolog as, without, almost any extensive Python package you want to use requires a Prolog library to make it all look acceptable. That means you still must write a comprehensive Prolog library and document it. My idea was to get the interface powerful enough to be able to make your application code directly use the Python package such that you no longer need a properly documented Prolog library but can use the Python docs right away.

Extensive use of e.g., numpy from Prolog probably still asks for a Prolog library. If you just need some isolated expensive operations on a couple of matrixes using numpy directly is good enough.

mgondan1 · December 17, 2023, 8:46pm

So, this would be the C++ binding to lapack (armadillo, to be precise):

#include <armadillo>
#include <SWI-cpp2.h>

using namespace arma;

struct Matrix;
static PL_blob_t matrix = PL_BLOB_DEFINITION(Matrix, "matrix");

struct Matrix : public PlBlob
{ mat m;

  explicit Matrix() : PlBlob(&matrix) { }

  explicit Matrix(const long rows, const long cols)
    : PlBlob(&matrix)
  { m.set_size(rows, cols);
    m.zeros();
  }

  explicit Matrix(mat am)
    : PlBlob(&matrix)
  { m = am;
  }

  PL_BLOB_SIZE

  bool write_fields(IOSTREAM *s, int flags) const override
  { for(uword i=0; i<m.n_rows; i++)
    { for(uword j=0; j<m.n_cols; j++)
        if(!Sfprintf(s, " %.3f", m(i, j)))
          return false;
      if(!Sfprintf(s, "\n"))
        return false;
    }
   if(!Sfprintf(s, "\n"))
      return false;
    return true;
  }
};

PREDICATE(matrix, 3)
{ auto ref = std::unique_ptr<PlBlob>(new Matrix(A1.as_int32_t(), A2.as_int32_t()));
  return A3.unify_blob(&ref);
}

PREDICATE(matmult, 3)
{ auto m1 = PlBlobV<Matrix>::cast_ex(A1, matrix);
  auto m2 = PlBlobV<Matrix>::cast_ex(A2, matrix);
  mat m = m1->m * m2->m;
  auto ref = std::unique_ptr<PlBlob>(new Matrix(m));
  return A3.unify_blob(&ref);
}

Compile with swipl-ld -shared plblas.cpp -o plblas.so -larmadillo

Test with

matthias@rechner:~/plblas$ swipl
Welcome to SWI-Prolog (threaded, 64 bits, version 9.0.4)
?- load_foreign_library(plblas).
true.

?- matrix(3, 2, A), matrix(2, 3, B), matmult(A, B, C).
A = <matrix>(0x55f4fc13fe40 0,000 0,000
 0,000 0,000
 0,000 0,000
),
B = <matrix>(0x55f4fc13ff30 0,000 0,000 0,000
 0,000 0,000 0,000
),
C = <matrix>(0x55f4fc140040 0,000 0,000 0,000
 0,000 0,000 0,000
 0,000 0,000 0,000
).

For now, the matrices are all zeros, but I just want to see how much typing is needed for a direct lapack interface without python support.

Am I doing something obviously wrong/unnecessary?

peter.ludemann · December 18, 2023, 1:48am

I’ve improved my sample code: DOC: Improve blob sample code by kamahen · Pull Request #65 · SWI-Prolog/packages-cpp · GitHub

And a few suggestions (I haven’t tested with armadillo because I haven’t installed all the prerequisites, so there might be some typos here):

I couldn’t find the type mat anywhere – is this Mat?
You might prefer this definition for the constructor that takes a mat as a parameter:

explicit Matrix(const mat& am) : PlBlob(&matrix), m(am) { }

Also, the constructor that uses rows and cols should probably be this (assuming it uses the Mat constructor); and there’s no need for const with uword, although it doesn’t hurt:

explict Matrix(uword rows, uword cols) ...

The write_fields() method is usually allowed to default, or only prints a bit more than the pointer info (e.g., you might want to print out the rows and cols information); if you want to pretty-print the output, you can define a predicate, e.g. matrix_portray(Stream, Blob), which you then use with something like:

multifile user:portray/1.
user:portray(Matrix) :- blob(Matrix, matrix), !, matrix_portray(current_output, Matrix).

PREDICATE(matrix_portray, 2)
{ auto ref = PlBlobV<MyBlob>::cast_ex(A2, my_blob);
  PlStream strm(A1, 0);
  ref->portray(strm);
  return true;
}

void Matrix::portray(PlStream& strm) const
{ strm.printf("Matrix(rows=%ul cols=%ul)", m.n_rows, m.n_cols);
  for(uword i=0; i<m.n_rows; i++)
    { for(uword j=0; j<m.n_cols; j++)
        strm.printf(s, " %.3f", m(i, j));
      strm.printf("\n");
   strm.printf("\n");
}

Note that instead of using Sfprintf(), you can use PlStream::printf() and there is no need to check the return code from the calls to printf().

mgondan1 · December 18, 2023, 9:38am

All right, thanks!

matthias@rechner:~/plblas$ swipl
Welcome to SWI-Prolog (threaded, 64 bits, version 9.0.4)
?- use_module(plblas).
true.
?- matrix(3, 2, _A), matrix(2, 3, _B), matmult(_A, _B, C).
C = Matrix(rows=3l cols=3l)
 0,000 0,000 0,000
 0,000 0,000 0,000
 0,000 0,000 0,000

This is the pl-file:

:- module(matrix, []).
:- load_foreign_library(plblas).
:- multifile user:portray/1.

user:portray(Matrix) :-
    blob(Matrix, matrix),
    !, matrix_portray(current_output, Matrix).

This is the cpp:

#include <armadillo>
#include <SWI-cpp2.h>

using namespace arma;
struct Matrix;

static PL_blob_t matrix = PL_BLOB_DEFINITION(Matrix, "matrix");

struct Matrix : public PlBlob
{ mat m;

  explicit Matrix()
    : PlBlob(&matrix) { }

  explicit Matrix(uword rows, uword cols)
    : PlBlob(&matrix)
  { m.set_size(rows, cols);
    m.zeros();
  }

  explicit Matrix(const mat& am)
    : PlBlob(&matrix)
  { m = am;
  }

  PL_BLOB_SIZE

  void portray(PlStream& strm) const;
};

PREDICATE(matrix, 3)
{ auto ref = std::unique_ptr<PlBlob>(new Matrix(A1.as_int32_t(), A2.as_int32_t()));
  return A3.unify_blob(&ref);
}

PREDICATE(matrix_portray, 2)
{ auto ref = PlBlobV<Matrix>::cast_ex(A2, matrix);
  PlTerm a1(A1);
  PlStream s(a1, SIO_OUTPUT);
  ref->portray(s);
  return true;
}

PREDICATE(matmult, 3)
{ auto m1 = PlBlobV<Matrix>::cast_ex(A1, matrix);
  auto m2 = PlBlobV<Matrix>::cast_ex(A2, matrix);
  mat m = m1->m * m2->m;
  auto ref = std::unique_ptr<PlBlob>(new Matrix(m));
  return A3.unify_blob(&ref);
}

void Matrix::portray(PlStream& s) const
{ s.printf("Matrix(rows=%ul cols=%ul)\n", m.n_rows, m.n_cols);
  for(uword i=0; i<m.n_rows; i++)
  {
    for(uword j=0; j<m.n_cols; j++)
      s.printf(" %.3f", m(i, j));

    s.printf("\n");
  }
}

I’ll setup a github repo for a regular swipl “pack”, in case anyone wants to join. If I am doing it alone, it will take a bit of time

PS. armadillo’s matrix type seems to be lowercase mat (?), see here Armadillo: C++ library for linear algebra & scientific computing - API Documentation

Boris · December 18, 2023, 11:49am

Hi, sorry for the silly question: how did you chose this library in particular? I am asking because I have previously needed some of the functionality and I always got lost in the multitude of options I had for “a C++ library for linear algebra”.

mgondan1 · December 18, 2023, 12:41pm

Correct if I am wrong, but isn’t Lapack de facto standard for linear algebra since decades?

It’s written in Fortran, though. And, I mean, most people in this forum obviously have big sympathy for dinosaur programming languages (e.g. Colmerauer 1972), but enough is enough. So the only thing needed is a C++ wrapper, and there you can basically choose between lapack++ (irregularly maintained) and armadillo.

Boris · December 18, 2023, 1:32pm

Yes, but as you say, LAPACK is not even C++. There are many other C++ libraries for linear algebra, with different trade-offs. So you narrowed it down by saying “C++ wrapper for LAPACK”; thank you, this is exactly the answer I was looking for.

mgondan1 · December 19, 2023, 9:21pm

I guess we can’t really support C++ templates, or can we?

peter.ludemann · December 19, 2023, 11:39pm

Do you mean that you want to do something like the following?

template <typename T>
struct Matrix : public PlBlob {
  mat<T> m;
  explicit Matrix<T>(const mat<T>& am)
    : PlBlob(&matrix), m(am) { }
};

It should be possible to do this, although you’d need to do something with the declaration for matrix (perhaps something with the preprocessor’s # and ## operators?):

static PL_blob_t matrix = PL_BLOB_DEFINITION(Matrix<concrete_type>, "matrix<concrete_type>");

Of course, you don’t have to use PL_BLOB_DEFINITION – after all, it merely generates a PL_blob_t initializer, so you could define something different … you might be able to reuse or extend class PlBlobV, which is used as a C wrapper for the C++ methods. If you can figure out a general way of doing this, we could add it to SWI-cpp2.h.

mgondan1 · December 21, 2023, 4:13pm

I see, thanks.

The thing is: For C++ templates, we need types, which Prolog hasn’t. Therefore, at some points, the different matrix types would need to be dispatched. Whether that happens “outside” or “inside” the blobs is probably not much of a difference. So, for now, I’ll move on with the cpp2 interface as it is.

peter.ludemann · December 21, 2023, 4:59pm

Templates usually use types, but don’t have to – they can use string or integers, for example. They’re also Turing complete (although actual implementations often have problems with complex templates).
Here’s an example of computing fibonacci numbers using templates: fibonacci template C++ · GitHub

mgondan1 · December 21, 2023, 5:19pm

Unification just like in Prolog, very nice. I wasn’t aware of this.

peter.ludemann · December 21, 2023, 5:51pm

Yes.
Everyone who uses C++ templates or Makefiles is using something that looks like a variant of Prolog.
My version of Greenspun’s Tenth Rule:

Any sufficiently complicated Lisp program contains an ad hoc, informally-specified, bug-ridden, slow implementation of half of Prolog. (Also applies to languages other than Lisp, but I’ve observed it a number of times with Lisp programs and textbooks.) (And also applies to Prolog - see Meta-Circular Interpreter.)

mgondan1 · December 29, 2023, 8:09pm

Ok. At some point, I need to dispatch the different types of blobs (say, matrices of integers and floats). Currently, I do this in prolog,

user:portray(Matrix) :- 
    blob(Matrix, matrix), 
    !, matrix_portray(current_output, Matrix).

user:portray(Column) :-
    blob(Column, column),
    !, column_portray(current_output, Column).

The example distinguishes matrices from column vectors, which may not even be necessary, but illustrates my point.

If I move this to C++, some kind of templates could be used. I guess I need your help, though, since I don’t fully understand PlBlobV

peter.ludemann · December 30, 2023, 3:35am

PlBlobV<blob_class> can be thought of as a way of generating struct PL_blob_t for your blob – the C++ template facility is used more-or-less as a macro and the methods are used by PL_BLOB_DEFINITION.

Your problem, if I understand correctly, is that you’re getting a term that’s an unknown blob type and want to process it. You can check the type using PlTerm::is_blob(), something like the following (untested code), which fails if the argument isn’t a blob of type “matrix” or “column”:

static PL_blob_t matrix = ...
static PL_blob_t column = ...

PREDICATE(my_portray, 1)
{ PL_blob_t *blob_type; 
  PlCheckFail(A1.is_blob(&blob_type));
  if ( blob_type == matrix )
  { ... matrix_portray ... 
  } else if ( blob_type == column )
  { ... column_portray ...
  } else
  { throw PlFail();
  }
  return true;
}

Boris · January 3, 2024, 8:44pm

Do you think it would be possible, using this library/binding, to define your own types and operations on them and apply the matrix or vector operations provided by armadillo? Like, could I give the dot product algorithm two vectors of some Prolog terms, and a “addition” and “multiplication” operation? Not so exciting for dot product but starts to get interesting for matrix multiplication.

EDIT: OK it sounds like a no I found this in the docs:

Is it possible to use Armadillo matrices with user-defined/custom element types ?
Armadillo supports matrices with the following element types: float, double, std::complex, std::complex, short, int, long, and unsigned versions of short, int, long. Support for other types is beyond the scope of Armadillo.

https://arma.sourceforge.net/faq.html#features

So I guess the answer is no unless I misunderstood something.

mgondan1 · January 5, 2024, 11:52am

Fine, thanks! Matrix and Column are now fully dispatched in C++. That is the relevant code:

static foreign_t plblas_portray(term_t a1, int arity, control_t ctx)
{ PL_blob_t* t ;
  if(((PlTerm) a1).is_blob(&t) == false)
    return false ;

  if(t == &matrix)
  { auto ref = PlBlobV<Matrix>::cast_ex((PlTerm) a1, matrix) ;
    PlStream s(PlTerm_atom("current_output"), SIO_OUTPUT) ;
    ref->portray(s) ;
    return true ;
  }

  if(t == &column)
  { auto ref = PlBlobV<Column>::cast_ex((PlTerm) a1, column) ;
    PlStream s(PlTerm_atom("current_output"), SIO_OUTPUT) ;
    ref->portray(s) ;
    return true ;
  }

  return false ;
}

PlRegister x_plblas_portray_1("user", "portray", 1, plblas_portray) ;

The prolog file of the pack only includes a :- multifile user:portray/1. I guess a few things can still be optimized, but none of them seem to be performance-critical.

The predicate has to return false to hand control back to prolog for portraying other objects than Matrix and Column blobs.

?- use_module(library(plblas)).
true.

?- matrix(2, 3, M).
M = Matrix(rows=2 cols=3)
 0.000 0.000 0.000
 0.000 0.000 0.000
.
?- column(3, M).
M = Column(rows=3)
 0.000
 0.000
 0.000
.
?- X = t(a).
X = t(a).

I’ll write back when I have cleaned up the other “blob methods”, and when I need help with the templates.

mgondan1 · January 8, 2024, 8:03am

I wasn’t aware that portray/1 is actually used in the backtrace:

?- use_module(library(plblas)).
true.

?- matrix(2, 2, A), eye(A).
A = Matrix(rows=2 cols=2)
 1.000 0.000
 0.000 1.000

.

?- column(2, A), eye(A).
ERROR: Type error: `Matrix' expected, found `Column(rows=2)
 0.000
 0.000

' (a column_reference)
ERROR: In:
ERROR:   [11] eye(Column(rows=2)
 0.000
 0.000

)
ERROR:    [9] toplevel_call(user:user: ...) at /usr/lib/swi-prolog/boot/toplevel.pl:1173
ERROR:
ERROR: Note: some frames are missing due to last-call optimization.
ERROR: Re-run your program in debug mode (:- debug.) to get more detail.

The response to matrix/eye is fine. The error message to column/eye is intended (because for now, eye/1 is reserved for matrices). It doesn’t look pretty though because the column is “portrayed” in the backtrace.

I guess I need to overwrite something else as portray (?)

jan · January 8, 2024, 9:05am

Toplevel answers as well as trace and backtrace output all use write_term/2 with the portray(true) option by default. There are several flags that allow controlling the options passed to write_term/2, e.g., answer_write_options, debugger_write_options and print_write_options. It is probably best to make sure your library works fine with the default settings.

I don’t think that printing a matrix content as default portray/1 behavior is a good idea. I think you want something terse for both the toplevel and debugger and a predicate to print the content of the matrix.

Following the normal blob conventions, the default print should be

 <matrix>(0xNNNN),

optionally adding the dimensions in some terse notation.

Topic		Replies	Views
Ann: SWI-Prolog 9.1.17 Releases	19	507	October 24, 2023
Swi-prolog bindings to SDL General	17	218	July 22, 2024
Install/use swi-prolog for C++ General	16	796	August 19, 2023
Ann: SWI-Prolog 9.0.0 Releases stable	1	1306	November 27, 2022
Changes to SWI-cpp.h Discussion discussion	62	2833	July 12, 2023

Lapack

Related topics