import scipy.io as sio

from scipy.sparse import lil_matrix
from scipy.sparse import identity

import numpy as np

#import sys
#if len(sys.argv) != 2:
#    print ("Usage: python ", sys.argv[0], " Calculixbasename")
#    quit()
#c_filename = sys.argv[1] + "_C.mtx"
#m_filename = sys.argv[1] + "_M.mtx"

m_filename = "FEMMeshNetgen" + "_M.mtx"
c_filename = "FEMMeshNetgen" + "_C.mtx"
dof_calculix_filname = "FEMMeshNetgen" + ".dof"


C_mat=lil_matrix(sio.mmread(c_filename))
M_mat=lil_matrix(sio.mmread(m_filename))


mat_dim = C_mat.get_shape()[1]
 
# #
# # check symmetry
# #
# v_one = np.ones(mat_dim)
# c_sym_err = (C_mat - C_mat.T) * (C_mat - C_mat.T) 
# print(c_sym_err.dot(v_one).dot(v_one))
# quit()


#D_mat = lil_matrix((mat_dim,mat_dim))
D_mat = C_mat * 0.1

import pandas as pd
dof_calculix = pd.read_csv(dof_calculix_filname,sep = '\.', engine='python',header=None, names=['Node', 'DOF'])


def calculix_dof2python_dof(dof_data, node, dof):
    query = '(Node == ' + str(node) + ' ) and (DOF == ' + str(dof) + ')'
    return(dof_data.query(query).index.item())

 
#
# vtk 943 - calculix 944
#

num_inputs = 1
B_mat = lil_matrix((mat_dim,num_inputs))
B_mat[calculix_dof2python_dof(dof_calculix,944,3) ,0] = 1.0

num_outputs = 1 
Cp = lil_matrix((num_outputs,mat_dim))
Cp[0,calculix_dof2python_dof(dof_calculix,3,3)] = 1.0


from pymor.models.iosys import SecondOrderModel
# Full Model
som = SecondOrderModel.from_matrices(M_mat, D_mat, C_mat, B_mat, Cp, Cv=None, D=None, 
    sampling_time=0, state_id='STATE', solver_options={ 'lyap_lrcf': {'type': 'lradi', 'tol': 1e-8}}, error_estimator=None, visualizer=None, name=None)


from  pymor.reductors.sobt import SOBTReductor
sobt = SOBTReductor(som)
rom = sobt.reduce(1)

rom.K.export_matrix('c_red.mtx',output_format='matrixmarket')
