1% ------------------------------------------------------
    2% Package: matrix_utls
    3% Title:   Matrix utilities package consists of:
    4%          1.  Kronecker (Tensor) product,
    5%          2.  Hadamard (Element-Wise) product,
    6%          3.  Creating (MxN) matrix of empty cells,
    7%          4.  Creating (MxN) matrix of a unique (int/float/letters) value,
    8%          5.  Creating (MxN) matrix of random (int/float) values,
    9%          6.  Creating (MxM) Identity matrix,
   10%          7.  Scalar-by-Vector multiplication,
   11%          8.  Scalar-by-Matrix multiplication,
   12%          9.  Vector-by-Vector multiplication,
   13%          10. Vector-by-Matrix multiplication,
   14%          11. Matrix-by-Vector multiplication, and
   15%          12. Matrix-by-Matrix multiplication.
   16% Version: 1.1
   17% Authors: Ali Al-Bayaty   <albayaty@pdx.edu>
   18%          Marek Perkowski <h8mp@pdx.edu>
   19%          Electrical & Computer Engineering Dept.
   20%          Portland State University
   21% Date:    02/25/2020
   22% ------------------------------------------------------
   23
   24
   25:- module(matrix_utls,
   26          [   kronecker/3,
   27              hadamard/3,
   28              create_empty_matrix/3,
   29              create_val_matrix/4,
   30              create_rand_matrix/5,
   31              create_I_matrix/2,
   32              s_v_multiply/3,
   33              s_m_multiply/3,
   34              v_v_multiply/3,
   35              v_m_multiply/3,
   36              m_v_multiply/3,
   37              m_m_multiply/3
   38          ]).
   39
   40
   41% Saving and restoring the temporary results (facts) to/from the KB:
   42:- dynamic outRes/1, subOutRes/1, addRes/1.
 kronecker(A, B, O)
Kronecker (Tensor) Product.

Where, A is the input (multiplicand) matrix of size [GxH], B is the input (multiplier) matrix of size [JxK], and O is the output matrix of size [(GxJ)x(HxK)].

   53kronecker(A, B, O):-
   54   assert(outRes([])),
   55   assert(subOutRes([])),
   56   maplist(get_A_Rows(B), A),
   57   retract(subOutRes(_)),
   58   retract(outRes(O)).
   59
   60% Getting the rows of A:
   61get_A_Rows(B, A_Row):-
   62   maplist(get_B_Rows(A_Row), B).
   63
   64% Getting the cells from 1 row of A:
   65get_A_Cells(B_Row, A_Cell):-
   66   maplist(get_B_Cells(A_Cell), B_Row).
   67
   68% Getting the rows of B:
   69get_B_Rows(A_Row, B_Row):-
   70   maplist(get_A_Cells(B_Row), A_Row),
   71   retract(outRes(Out)),
   72   retract(subOutRes(SubOut)),
   73   append(Out, [SubOut], FinalOut),
   74   assert(subOutRes([])),
   75   assert(outRes(FinalOut)).
   76
   77% Getting the cells from 1 row of B:
   78get_B_Cells(A_Cell, B_Cell):-
   79   % Element-wise cells multiplications:
   80   R is A_Cell * B_Cell,
   81   retract(subOutRes(SubOut)),
   82   append(SubOut, [R], FinalSubOut),
   83   assert(subOutRes(FinalSubOut)).
 hadamard(A, B, O)
Hadamard (Element-Wise) Product.

Where, A is the input (multiplicand) matrix of size [GxH], B is the input (multiplier) matrix of size [GxH], and O is the output matrix of size [GxH].

   94hadamard(A, B, O):-
   95   % Checking the dimensional equality:
   96   checkRows(A ,B, Num_Rows),
   97   checkCols(A, B, Num_Cols),
   98   assert(outRes([])),
   99   assert(subOutRes([])),
  100   forall(between(1, Num_Rows, Row_Index),
  101                                   doRows(A, B, Row_Index, Num_Cols)),
  102   retract(subOutRes(_)),
  103   retract(outRes(O)).
  104
  105% Checking the rows equality:
  106checkRows(A, B, Num_A_Rows):-
  107   length(A, Num_A_Rows),
  108   length(B, Num_B_Rows),
  109   ( Num_A_Rows =\= Num_B_Rows
  110      -> writeln('\nError: Matrices do not have equal dimensionality (rows) ...\n'),
  111         abort
  112      ;  ! ).
  113
  114% Checking the cols equality:
  115checkCols(A, B, Num_A_Cols):-
  116   nth1(1, A, A_Row),
  117   nth1(1, B, B_Row),
  118   length(A_Row, Num_A_Cols),
  119   length(B_Row, Num_B_Cols),
  120   ( Num_A_Cols =\= Num_B_Cols
  121      -> writeln('\nError: Matrices do not have equal dimensionality (columns) ...\n'),
  122         abort
  123      ;  !).
  124
  125doRows(A, B, Row_Index, Num_Cols):-
  126   nth1(Row_Index, A, A_Row),
  127   nth1(Row_Index, B, B_Row),
  128   forall(between(1, Num_Cols, Col_Index),
  129                                    doCols(A_Row, B_Row, Col_Index)),
  130   retract(outRes(Out)),
  131   retract(subOutRes(SubOut)),
  132   append(Out, [SubOut], FinalOut),
  133   assert(subOutRes([])),
  134   assert(outRes(FinalOut)).
  135
  136doCols(A_Row, B_Row, Col_Index):-
  137   nth1(Col_Index, A_Row, A_Cell),
  138   nth1(Col_Index, B_Row, B_Cell),
  139   R is A_Cell * B_Cell,
  140   retract(subOutRes(SubOut)),
  141   append(SubOut, [R], FinalSubOut),
  142   assert(subOutRes(FinalSubOut)).
 create_empty_matrix(Rows, Cols, Mat)
Create 2D matrix of empty cells.

Where, Rows is the number of rows, Cols is the number of columns, and Mat is the resultant matrix of size [Rows x Cols].

  153create_empty_matrix(Rows, Cols, Mat):-
  154   length(Mat, Rows),
  155   % Creating temporary Lists based on the number of rows:
  156   maplist(create_list(Cols), Mat).
 create_val_matrix(Rows, Cols, Mat, Value)
Create 2D matrix and set with a unique (int/float/letters) value.

Where, Rows is the number of rows, Cols is the number of columns, and Mat is the resultant matrix of size [Rows x Cols].

  167create_val_matrix(Rows, Cols, Mat, Value):-
  168   create_empty_matrix(Rows, Cols, Mat),
  169   setValue(Value, 1, Mat).
 create_rand_matrix(Rows, Cols, Mat, Min_Value, Max_Value)
Create 2D matrix and set with a range of random (int/float) values.

Where, Rows is the number of rows, Cols is the number of columns, and Mat is the resultant matrix of size [Rows x Cols].

  180create_rand_matrix(Rows, Cols, Mat, Min_Value, Max_Value):-
  181   create_empty_matrix(Rows, Cols, Mat),
  182   setValue([Min_Value, Max_Value], 2, Mat).
  183
  184% Creating the columns of Mat:
  185create_list(Cols, Sub_Mat):-
  186   length(Sub_Mat, Cols).
  187
  188% Setting values for each cell of Mat:
  189setValue(Value, Type, Mat):-
  190   maplist(setRows(Value, Type), Mat).
  191setRows(Value, Type, Mat_Rows):-
  192   maplist(setCells(Value, Type), Mat_Rows).
  193setCells(Value, 1, Mat_Cell):- Mat_Cell = Value, !.
  194setCells(Value, 2, Mat_Cell):-
  195   nth1(1, Value, Min_Value),
  196   nth1(2, Value, Max_Value),
  197   random(Min_Value, Max_Value, Mat_Cell), !.
 create_I_matrix(Dim, Mat)
Creating Identity (MxM) Matrix.

Where, Dim is the number of rows and columns, and Mat is the resultant identity matrix of size [Dim x Dim].

  207create_I_matrix(Dim, Mat):-
  208   Dim_1 is Dim - 1,
  209   create_val_matrix(1, Dim_1, Zeros, 0),
  210   nth1(1, Zeros, Zero),
  211   length(Mat, Dim),
  212   numlist(1, Dim, Index),
  213   maplist(insert_Ones(Zero), Index, Mat).
  214
  215insert_Ones(Zero, Index, Sub_Mat):-
  216   nth1(Index, Sub_Mat, 1, Zero).
 s_v_multiply(S, V, O)
Scalar-by-Vector Multiplication.

Where, S is the input (multiplicand) scalar (i.e. number), V is the input (multiplier) vector of size [1xH] or [Gx1], and O is the resultant vector of size [1xH] or [Gx1].

  227s_v_multiply(S, V, O):-
  228   s_m_multiply(S, V, O).
 s_m_multiply(S, M, O)
Scalar-by-Matrix Multiplication.

Where, S is the input (multiplicand) scalar (i.e. number), M is the input (multiplier) matrix of size [GxH], and O is the resultant matrix of size [GxH].

  239s_m_multiply(S, M, O):-
  240   assert(outRes([])),
  241   assert(subOutRes([])),
  242   length(M, M_Rows),
  243   forall(between(1, M_Rows, Row_Index), do_Rows(S, M, Row_Index)),
  244   retract(subOutRes(_)),
  245   retract(outRes(O)).
  246
  247do_Rows(S, M, Row_Index):-
  248   nth1(Row_Index, M, M_Row),
  249   length(M_Row, M_Cols),
  250   forall(between(1, M_Cols, Col_Index), do_Cols(S, M_Row, Col_Index)),
  251   retract(outRes(Out)),
  252   retract(subOutRes(SubOut)),
  253   append(Out, [SubOut], FinalOut),
  254   assert(subOutRes([])),
  255   assert(outRes(FinalOut)).
  256
  257do_Cols(S, M_Row, Col_Index):-
  258   nth1(Col_Index, M_Row, M_Cell),
  259   R is S * M_Cell,
  260   retract(subOutRes(SubOut)),
  261   append(SubOut, [R], FinalSubOut),
  262   assert(subOutRes(FinalSubOut)).
 v_v_multiply(V1, V2, O)
Vector-by-Vector Multiplication.

Where, V1 is the input (multiplicand) vector of size [1xJ] or [Kx1], V2 is the input (multiplier) vector of size [Jx1] or [1xK], and O is the resultant scalar of size [1x1] or matrix of size [KxK].

  273v_v_multiply(V1, V2, O):-
  274   m_m_multiply(V1, V2, O).
 v_m_multiply(V, M, O)
Vector-by-Matrix Multiplication.

Where, V is the input (multiplicand) vector of size [1xJ], M is the input (multiplier) matrix of size [JxK], and O is the resultant vector of size [1xK].

  285v_m_multiply(V, M, O):-
  286   m_m_multiply(V, M, O).
 m_v_multiply(M, V, O)
Matrix-by-Vector Multiplication.

Where, M is the input (multiplicand) matrix of size [JxK], V is the input (multiplier) vector of size [Kx1], and O is the resultant vector of size [Jx1].

  297m_v_multiply(M, V, O):-
  298   m_m_multiply(M, V, O).
 m_m_multiply(M1, M2, O)
Matrix-by-Matrix Multiplication.

Where, M1 is the input (multiplicand) matrix of size [GxH], M2 is the input (multiplier) matrix of size [HxJ], and O is the resultant matrix of size [GxJ].

  309m_m_multiply(M1, M2, O):-
  310   length(M1, Num_M1_Rows),
  311   length(M2, Num_M2_Rows),
  312   nth1(1, M1, M1_Row),
  313   length(M1_Row, Num_M1_Cols),
  314   nth1(1, M2, M2_Row),
  315   length(M2_Row, Num_M2_Cols),
  316   checkDim(Num_M1_Cols, Num_M2_Rows),
  317   assert(outRes([])),
  318   assert(subOutRes([])),
  319   assert(addRes(0.0)),
  320   forall(between(1, Num_M1_Rows, M1_Row_Index),
  321          do_M1_Rows(M1, M2, M1_Row_Index, Num_M1_Cols, Num_M2_Cols)),
  322   retract(addRes(_)),
  323   retract(subOutRes(_)),
  324   retract(outRes(O)).
  325
  326% Checking the equal dimensionality requirement:
  327checkDim(Num_M1_Cols, Num_M2_Rows):-
  328   ( Num_M1_Cols =\= Num_M2_Rows
  329      -> writeln('\nError: Matrices do not have equal dimensionality requirement:'),
  330         write('M1 cols = '), write(Num_M1_Cols), nl,
  331         write('M2 rows = '), write(Num_M2_Rows), nl,
  332         writeln('Hint: Both matrices should have equal dimensions (rows and cols).\n'),
  333         abort
  334      ;  ! ).
  335
  336do_M1_Rows(M1, M2, M1_Row_Index, Num_M1_Cols, Num_M2_Cols):-
  337   forall(between(1, Num_M2_Cols, M2_Col_Index),
  338          do_M2_Col(M1, M2, M1_Row_Index,Num_M1_Cols, M2_Col_Index)),
  339   retract(outRes(Out)),
  340   retract(subOutRes(SubOut)),
  341   append(Out, [SubOut], FinalOut),
  342   assert(subOutRes([])),
  343   assert(outRes(FinalOut)).
  344
  345do_M2_Col(M1, M2, M1_Row_Index,Num_M1_Cols, M2_Col_Index):-
  346   nth1(M1_Row_Index, M1, M1_Row),
  347   forall(between(1, Num_M1_Cols, M1_Col_Index),
  348          doMultiply(M1_Row, M2, M1_Col_Index, M2_Col_Index)),
  349   retract(subOutRes(SubOut)),
  350   retract(addRes(SubAdd)),
  351   append(SubOut, [SubAdd], FinalSubOut),
  352   assert(addRes(0.0)),
  353   assert(subOutRes(FinalSubOut)).
  354
  355doMultiply(M1_Row, M2, M1_Col_Index, M2_Col_Index):-
  356   nth1(M1_Col_Index, M1_Row, M1_Cell),
  357   nth1(M1_Col_Index, M2, M2_Row),
  358   nth1(M2_Col_Index, M2_Row, M2_Cell),
  359   R is M1_Cell * M2_Cell,
  360   retract(addRes(SubAdd)),
  361   Out is SubAdd + R,
  362   assert(addRes(Out)).
  363
  364
  365% ------------------------------------------------------
  366% %%%%%%%%%%%%%%%%%%% EXAMPLES %%%%%%%%%%%%%%%%%%%%%%%%%
  367% ------------------------------------------------------
  368% ?- kronecker([[1,2],[2,-1]], [[1,2],[3,4]], M).
  369%    M = [[1, 2,  2,  4],
  370%         [3, 4,  6,  8],
  371%         [2, 4, -1, -2],
  372%         [6, 8, -3, -4]].
  373%
  374% ?- hadamard([[1,2],[2,-1]], [[1,2],[3,4]], M).
  375%    M = [[1,  4],
  376%         [6, -4]].
  377%
  378% ?- create_empty_matrix(2, 3, M).
  379%    M = [[_928, _934, _940],
  380%         [_946, _952, _958]].
  381%
  382% ?- create_val_matrix(2, 3, M, 1).
  383%    M = [[1, 1, 1],
  384%         [1, 1, 1]].
  385%
  386% ?- create_val_matrix(2, 3, M, Aa).
  387%    M = [[Aa, Aa, Aa],
  388%         [Aa, Aa, Aa]].
  389%
  390% ?- create_rand_matrix(2, 3, M, -2.0, 2.0).
  391%    M = [[-0.331722127869,  1.2472310874774, -0.3789740703109],
  392%         [0.04395141605513, 0.5567944218039, -0.27190713107775]].
  393%
  394% ?- create_I_matrix(3, M).
  395%    M = [[1, 0, 0],
  396%         [0, 1, 0],
  397%         [0, 0, 1]].
  398%
  399% ?- s_v_multiply(5, [[1, 2, 3]], M).
  400%    M = [[5, 10, 15]].
  401%
  402% ?- s_v_multiply(5, [[1], [2], [3]], M).
  403%    M = [[5],
  404%         [10],
  405%         [15]].
  406%
  407% ?- s_m_multiply(5, [[1,2,3], [4,5,6]], M).
  408%    M = [[5,  10, 15],
  409%         [20, 25, 30]].
  410%
  411% ?- v_v_multiply([[1, 2, 3]], [[1], [2], [3]], M).
  412%    M = [[14.0]].
  413%
  414% ?- v_v_multiply([[1], [2], [3]], [[1, 2, 3]], M).
  415%    M = [[1.0, 2.0, 3.0],
  416%         [2.0, 4.0, 6.0],
  417%         [3.0, 6.0, 9.0]].
  418%
  419% ?- v_m_multiply([[1, 2, 3]], [[1,2], [3,4], [5,6]], M).
  420%    M = [[22.0, 28.0]].
  421%
  422% ?- m_v_multiply([[1,2], [3,4], [5,6]], [[1], [2]], M).
  423%    M = [[5.0],
  424%         [11.0],
  425%         [17.0]].
  426%
  427% ?- m_m_multiply([[1,2,3],[4,5,6],[7,8,9]],
  428%                 [[1,0,0],[0,1,0],[0,0,1]], M).
  429%    M = [[1.0, 2.0, 3.0],
  430%         [4.0, 5.0, 6.0],
  431%         [7.0, 8.0, 9.0]].
  432%
  433% ------------------------------------------------------