1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
|
import itertools
import numpy as np
from polymatrix.statemonad.init import init_state_monad
from polymatrix.statemonad.mixins import StateMonadMixin
from polymatrix.expressionstate import ExpressionState
from polymatrix.expression.expression import Expression
from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
from polymatrix.expression.utils.getvariableindices import (
get_variable_indices_from_variable,
)
from polymatrix.denserepr.utils.monomialtoindex import monomial_to_monomial_vector_indices
from polymatrix.denserepr.impl import DenseReprBufferImpl, DenseReprImpl
# NP: create a dense representation from a polymatrix expression
# FIXME:
# - typing problems
# - create custom exception for error (instead of AssertionError & Exception)
def from_polymatrix(
expressions: Expression | tuple[Expression],
variables: Expression = None,
sorted: bool = None,
) -> StateMonadMixin[
ExpressionState, tuple[tuple[tuple[np.ndarray, ...], ...], tuple[int, ...]]
]:
"""
Converts a tuple of polynomial vectors `expressions` into a (semi-dense) matrix representation.
Given single polynomial vector
expr = [[x**2 + y], [1 + x*y + y**2]],
its matrix representation is given by three matrices
1. np.array type matrix [[0], [1]] containing the constant part
2. np.array type matrix [[0, 1], [0, 0]] containing the linear part
w.r.t the monomial vector Z=[x, y]
3. scipy.sparse type matrix [[1, 0, 0, 0], [0, 0.5, 0.5, 1]] containing the quadratic part
w.r.t the monomial vector Z=[x**2, x*y, x*y, y**2].
Instead of providing only a single polynomial vector, a tuple of polynomial vectors can be provided.
This approach ensures consistency in the vector of monomial accross all matrix representations.
"""
if isinstance(expressions, Expression):
expressions = (expressions,)
assert (
isinstance(variables, ExpressionBaseMixin) or variables is None
), f"{variables=}"
def func(state: ExpressionState):
def acc_underlying_application(acc, v):
state, underlying_list = acc
state, underlying = v.apply(state)
assert underlying.shape[1] == 1, f"{underlying.shape[1]=} is not 1"
return state, underlying_list + (underlying,)
*_, (state, polymatrix_list) = tuple(
itertools.accumulate(
expressions,
acc_underlying_application,
initial=(state, tuple()),
)
)
if variables is None:
sorted_variable_index = tuple()
else:
state, variable_index = get_variable_indices_from_variable(state, variables)
if sorted:
tagged_variable_index = tuple(
(offset, state.get_name_from_offset(offset))
for offset in variable_index
)
sorted_variable_index = tuple(
v[0] for v in sorted(tagged_variable_index, key=lambda v: v[1])
)
else:
sorted_variable_index = variable_index
sorted_variable_index_set = set(sorted_variable_index)
if len(sorted_variable_index) != len(sorted_variable_index_set):
duplicates = tuple(
state.get_name_from_offset(var)
for var in sorted_variable_index_set
if 1 < sorted_variable_index.count(var)
)
raise Exception(
f"{duplicates=}. Make sure you give a unique name for each variables."
)
# This dictionary maps the expression variable index (saved in the expression state)
# to the variable index of the dense representation.
# Assume state.indices = {'x': (2,3), 'y': (1,2), 'z': (0,1)} and variable = [[2], [1]],
# then variable_index_map = {2: 0, 1: 1}.
variable_index_map = {old: new for new, old in enumerate(sorted_variable_index)}
n_param = len(sorted_variable_index)
def gen_numpy_matrices():
for polymatrix in polymatrix_list:
n_row = polymatrix.shape[0]
buffer = DenseReprBufferImpl(
data={},
n_row=n_row,
n_param=n_param,
)
# MS: ignores any columns, better would be to check that there are no columns
for row in range(n_row):
polymatrix_terms = polymatrix.get_poly(row, 0)
if polymatrix_terms is None:
continue
if len(polymatrix_terms) == 0:
buffer.add(row, 0, 0, 0)
else:
for monomial, value in polymatrix_terms.items():
def gen_new_monomial():
for var, count in monomial:
try:
index = variable_index_map[var]
except KeyError:
# todo: improve this error message!
raise KeyError(
f"{var=} ({state.get_key_from_offset(var)}) is incompatible with {variable_index_map=}"
)
for _ in range(count):
yield index
# converts the monomial x**2 = ((2, 2),) to (2, 2)
# converts the monomial x*y = ((2, 1), (1, 1)) to (2, 1)
# converts the monomial y**2 = ((1, 2),) to (1, 1)
new_variable_indices = tuple(gen_new_monomial())
# converts (2, 2) to (0,)
# converts (2, 1) to (1, 2)
# converts (1, 2) to (3,)
# the cols correspond to the column of the dense matrix w.r.t. Z
cols = monomial_to_monomial_vector_indices(
n_param, new_variable_indices
)
# the monomial x*y is mapped to two columns, hence we divide its value by 2:
# x*y = [0, 0.5, 0.5, 0] @ [x**2, x*y, x*y, y**2].T
col_value = value / len(cols)
for col in cols:
degree = sum(count for _, count in monomial)
buffer.add(row, col, degree, col_value)
yield buffer
underlying_matrices = tuple(gen_numpy_matrices())
result = DenseReprImpl(
data=underlying_matrices,
variable_mapping=sorted_variable_index,
state=state,
)
return state, result
return init_state_monad(func)
|
