summaryrefslogtreecommitdiffstats
path: root/polymatrix/denserepr/from_.py
blob: 13b77b76ae267bcd1496d0f5c663ddf20670beee (plain)
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
import dataclasses
import itertools
import typing
import numpy as np
import scipy.sparse
from polymatrix.denserepr.impl import DenseReprBufferImpl, DenseReprImpl

from polymatrix.expression.expression import Expression
from polymatrix.expression.mixins.expressionbasemixin import ExpressionBaseMixin
from polymatrix.expressionstate.abc import ExpressionState
from polymatrix.expression.utils.getvariableindices import get_variable_indices_from_variable
from polymatrix.statemonad.init.initstatemonad import init_state_monad
from polymatrix.statemonad.mixins.statemonadmixin import StateMonadMixin
from polymatrix.expression.utils.monomialtoindex import monomial_to_index


def from_polymatrix_expr(
    expressions: Expression | tuple[Expression],
    variables: Expression = None,
    sorted: bool = None,
) -> StateMonadMixin[ExpressionState, tuple[tuple[tuple[np.ndarray, ...], ...], tuple[int, ...]]]:

    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.')
        
        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,
                )

                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:
                                        raise KeyError(f'{var=} ({state.get_key_from_offset(var)}) is incompatible with {variable_index_map=}')

                                    for _ in range(count):
                                        yield index

                            new_monomial = tuple(gen_new_monomial())

                            cols = monomial_to_index(n_param, new_monomial)

                            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())

        def gen_auxillary_equations():
            for key, monomial_terms in state.auxillary_equations.items():
                if key in sorted_variable_index:
                    yield key, monomial_terms

        auxillary_equations = tuple(gen_auxillary_equations())

        n_row = len(auxillary_equations)

        if n_row == 0:
            auxillary_matrix_equations = None

        else:
            buffer = DenseReprBufferImpl(
                data={},
                n_row=n_row,
                n_param=n_param,
            )

            for row, (key, monomial_terms) in enumerate(auxillary_equations):
                for monomial, value in monomial_terms.items():
                    new_monomial = tuple(variable_index_map[var] for var, count in monomial for _ in range(count))

                    cols = monomial_to_index(n_param, new_monomial)

                    col_value = value / len(cols)

                    for col in cols:
                        buffer.add(row, col, sum(count for _, count in monomial), col_value)

            auxillary_matrix_equations = buffer

        result = DenseReprImpl(
            data=underlying_matrices,
            aux_data=auxillary_matrix_equations,
            variable_mapping=sorted_variable_index,
            state=state,
        )

        return state, result

    return init_state_monad(func)