summaryrefslogtreecommitdiffstats
path: root/main.py
diff options
context:
space:
mode:
authorMichael Schneeberger <michael.schneeberger@fhnw.ch>2022-04-12 09:30:35 +0200
committerMichael Schneeberger <michael.schneeberger@fhnw.ch>2022-04-12 09:30:35 +0200
commitf9ec423683462800c757119c08947b742458639e (patch)
treec75a4686696368b160a46ed4f67601be8581492f /main.py
parentremove shape property, introduce accumulate function, reimplement derivative ... (diff)
downloadpolymatrix-f9ec423683462800c757119c08947b742458639e.tar.gz
polymatrix-f9ec423683462800c757119c08947b742458639e.zip
bug fixes and clean ups
Diffstat (limited to 'main.py')
-rw-r--r--main.py159
1 files changed, 0 insertions, 159 deletions
diff --git a/main.py b/main.py
deleted file mode 100644
index 4079ac5..0000000
--- a/main.py
+++ /dev/null
@@ -1,159 +0,0 @@
-import sympy
-import scipy.special
-import scipy.sparse
-import numpy as np
-from polymatrix.polystruct import init_equation, init_poly_matrix, init_poly_vector
-from polymatrix.sympyutils import poly_to_data_coord, poly_to_matrix
-from polymatrix.utils import variable_to_index
-
-binom = scipy.special.binom
-
-def skew_symmetric(degree, p_row, p_col):
- if p_col < p_row:
- return p_col, p_row, -1
-
-def zero_diagonal_for_uneven_degree(degree, p_row, p_col):
- if degree % 2 == 1 and p_row == p_col:
- return 0
-
-def gradient(degree, v_row, monom):
- if degree == 1:
- factor = sum(v_row==e for e in monom) + 1
-
- if monom[-1] < v_row:
- n_v_row = monom[-1]
- n_monom = sorted(monom + (v_row,), reverse=True)
-
- if v_row <= monom[-1]:
- n_v_row = v_row
- n_monom = monom
-
- return n_v_row, n_monom, factor
-
-x = sympy.symbols('v_dc, i_d, i_q')
-vdc, id, iq = x
-
-g_dc = 3
-wl = 106
-
-jr_poly_list = [
- g_dc, -1, g_dc*wl,
- 1, 0, wl,
- -g_dc*wl, -wl, 0
-]
-jr_terms = poly_to_data_coord(list(sympy.poly(p, x) for p in jr_poly_list))
-# print(jr_terms)
-
-h = vdc**2/2 + id**2/2 + iq**2/2
-dh_poly_list = [
- sympy.diff(h, vdc), sympy.diff(h, id), sympy.diff(h, iq)
-]
-dh_terms = poly_to_data_coord(list(sympy.poly(p, x) for p in dh_poly_list))
-# print(dh_terms)
-
-mg_poly_list = [
- g_dc-id, -iq,
- vdc+1, 0,
- 0, vdc+1
-]
-mg_terms = poly_to_data_coord(list(sympy.poly(p, x) for p in mg_poly_list))
-# print(mg_terms)
-
-nabla_h = init_poly_vector(subs=dh_terms)
-JR = init_poly_matrix(subs=jr_terms)
-mg = init_poly_matrix(subs=mg_terms)
-
-nabla_ha = init_poly_vector(
- degrees=(1,),
- re_index_func_2=gradient,
-)
-JRa = init_poly_matrix(
- degrees=(0,1,2),
- re_index_func=skew_symmetric,
- subs_func=zero_diagonal_for_uneven_degree,
-)
-u = init_poly_vector(degrees=(1,2))
-
-eq = init_equation(
- terms = [(JR, nabla_ha), (JRa, nabla_ha), (JRa, nabla_h), (mg, u)],
- n_var = 2,
- # terms = [(JR, nabla_ha)],
- # n_var = 2,
-)
-
-# n_var = 2
-# n_param_nabla_ha = sum(n_var*binom(n_var+p-1, p) for p in nabla_ha.degrees)
-# n_param_JRa = sum(n_var**2*binom(n_var+p-1, p) for p in JRa.non_zero_degrees)
-# n_param_u = sum(n_var*binom(n_var+p-1, p) for p in u.non_zero_degrees)
-# total = n_param_nabla_ha+n_param_JRa+n_param_u
-
-# print(f'{n_param_nabla_ha=}')
-# print(f'{n_param_JRa=}')
-# print(f'{n_param_u=}')
-# print(f'{total=}')
-
-# print(binom(total+2-1, 2))
-
-# mat = init_poly_matrix(
-# degrees=(1,),
-# re_index_func=skew_symmetric,
-# subs_func=zero_diagonal_for_uneven_degree,
-# )
-# vec = init_poly_vector(
-# subs={0: {(0, 0): 1, (1, 0): 1}},
-# )
-
-# # mat = init_poly_matrix(
-# # subs={0: {(0, 0): 1, (1, 0): 1, (0, 1): 1, (1, 1): 1}},
-# # )
-
-# # vec = init_poly_vector(
-# # degrees=(1,),
-# # re_index_func_2=gradient,
-# # )
-
-# eq = init_equation(
-# terms = [(mat, vec)],
-# n_var = 2,
-# )
-
-eq_tuples, offset_dict, n_param = eq.create()
-
-print(f'{list(offset_dict.values())=}')
-print(f'{n_param=}')
-
-# mapping from row index to entry in eq_tuples
-rows_to_eq = list(set(key for eq_tuple_degree in eq_tuples.values() for key in eq_tuple_degree.keys()))
-eq_to_rows = {eq: idx for idx, eq in enumerate(rows_to_eq)}
-
-print(f'n_eq = {len(rows_to_eq)}')
-
-# # mapping from col index to entry in eq_tuples
-# cols_to_var = list(set(key for eq in eq_tuples[degree].values() for key in eq.keys()))
-# var_to_cols = {var: idx for idx, var in enumerate(cols_to_var)}
-
-def gen_matrix_per_degree():
- for degree, eq_tuple_degree in eq_tuples.items():
-
- if 0 < degree:
- def gen_coords(eq_tuple_degree=eq_tuple_degree):
- for eq, var_dict in eq_tuple_degree.items():
- zero_degree_val = eq_tuples[0][eq][0]
- row = eq_to_rows[eq]
- for var, value in var_dict.items():
- # col = variable_to_index(n_param, var)
- # col = var_to_cols[var]
- col = var
- yield row, col, value / zero_degree_val
-
- row, col, data = list(zip(*gen_coords()))
-
- n_col = int(binom(n_param+degree-1, degree))
- n_row = int(max(row) + 1)
- sparse_matrix = scipy.sparse.coo_matrix((data, (row, col)), dtype=np.float64, shape=(n_row, n_col)).toarray()
-
- yield degree, sparse_matrix
-
-result = dict(gen_matrix_per_degree())
-
-print(result[1].shape)