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

 

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-3}}, error_estimator=None, visualizer=None, name=None)


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

#from  pymor.reductors.sobt import SOBTfvReductor
#sobt = SOBTfvReductor(som)
#rom = sobt.reduce(20)
#
# same as SOBTReductor
#
# ------------------------------------------------------------------------


#from pymor.reductors.sor_irka import SORIRKAReductor
#sobt = SORIRKAReductor(som)
#rom = sobt.reduce(20)
#
# 00:00 SORIRKAReductor: Generating initial interpolation data
# 00:00 SORIRKAReductor: Starting SOR-IRKA
#
# interrupted after some time

# ------------------------------------------------------------------------


# no change .....
#rom = sobt.reduce(20,projection='orth',tol=1e-2,num_prev=5)

# ------------------------------------------------------------------------


#lti=som.to_lti()
#from pymor.reductors.bt import BTReductor
#bt=BTReductor(lti)
#rom = bt.reduce(20)
#
#
# same as SOBTReductor, which makes sense
#

# ------------------------------------------------------------------------


#lti=som.to_lti()
#from pymor.reductors.bt import FDBTReductor
#bt=FDBTReductor(lti)
#rom = bt.reduce(20)
#
#00:00 |WARNING|LTIModel: Converting operator A to a NumPy array.
#00:00 |WARNING|LTIModel: Converting operator E to a NumPy array.
#
# converting sparse to dense matrices, does not work for FEA at reasonable number of DOF 
#

# ------------------------------------------------------------------------


#lti=som.to_lti()
#from pymor.reductors.bt import LQGBTReductor
#bt=LQGBTReductor(lti)
#rom = bt.reduce(20)
#
#
# same as SOBTReductor
#



print(rom)



quit()
