(* 
	Evgenii Rudnyi (c) 2004, 2005, 2006, http://Evgenii.Rudnyi.Ru/

	Functions to research on model reduction before implementing them in mor4ansys
	
*)

BeginPackage["ModelReduction`", {"Imtek`Post4MOR`", "Utilities`FilterOptions`"}];

ArnoldiProcess::usage = "ArnoldiProcess[nextVec, startVec, dim] gives transposed projection matrix V with dim columns generated by the Arnoldi process with the starting vector startVec and the function nextVec to generate the next vector as a function of a previous vector.";

ProjectSystem::usage = "ProjectSystem[sys, V] gives a reduced dynamic system projected on the basis V. ProjectSystem[psys, V] gives a reduced parametric system. V must be the transposed projection matrix as returned by ArnoldiProcess.";

ModelReduction::usage = "ModelReduction[sys, dim, ops] gives a reduced dynamic system with dimension dim. Options are ExpansionPoint and options for LinearSolve.";

ExpansionPoint::usage = "An option for ModelReduction (default is 0).";

Orthogonalize::usage = "Orthogonalize[Vini, V, tol] gives a joint orthogonalized basis of Vini and V. tol is the tolerance to deflate new vectors (by default 1*^-12).";

ParametricSystem::usage = "An object that represent a system for parametric model reduction. Use MakeParametricSystem to create it and Matrices, Matrix[i], MatrixB, MatrixC, OutputNames, TransformFunctions, MatrixMFunction, MatrixEFunction, MatrixKFunction to access individual members.";

MakeParametricSystem::usage = "MakeParametricSystem[mats, matB, matC, outNames, funs] gives a ParametricSystem. mats is a list of matrices and funs is a list of three functions that transforms mats to matM, matE and matK of DynamicSystem. For example, matM = funs[[1]][mats, par] and so on. Use Null& if you do not need a particular matrix of DynamicSystem.";

Matrices::usage = "Matrices[psys] gives system matrices of psys.";

Matrix::usage = "Matrix[psys, i] gives the i-th system matrix of psys.";

MatrixB::usage = "MatrixB[psys] gives matB of psys.";

MatrixC::usage = "MatrixC[psys] gives matC of psys.";

OutputNames::usage = "OutputNames[psys] gives outNames of psys.";

TransformFunctions::usage = "TransformFunctions[psys] gives funs of psys.";

MatrixMFunction::usage = "MatrixMFunction[psys] gives the tranformation matrix to obtain matM = MatrixMFunction[psys][mats, par].";

MatrixEFunction::usage = "MatrixEFunction[psys] gives the tranformation matrix to obtain matE = MatrixEFunction[psys][mats, par].";

MatrixKFunction::usage = "MatrixKFunction[psys] gives the tranformation matrix to obtain matK = MatrixKFunction[psys][mats, par].";

MakeDynamicSystem::usage = "MakeDynamicSystem[psys, {par}] gives DynamicSystem obtained from psys at values of parameters specified by a list.";

(*::usage = "";*)

Begin["`Private`"];

ArnoldiProcess[nextVector_Function, vec_List, size_Integer] := 
Module[{matrixV, vectorW, norm},
	matrixV = {vec/Norm[vec]};
	For[j = 1, j < size, ++j,
		vectorW = nextVector[matrixV[[j]]];
		For[i = 1, i <= j, ++i,
			vectorW = vectorW - (Conjugate[matrixV[[i]]].vectorW)*matrixV[[i]]
		];
		If[(norm = Norm[vectorW]) > 1*^-12,
		AppendTo[matrixV, vectorW/norm],
		Print["Deflation during vector ", j + 1]; Break[]]
	];
	matrixV
]

projectMatrix[mat_, V_List] := Module[{tmp},
	tmp = mat.Transpose[V];
	Conjugate[V].tmp
]

ProjectSystem[sys_DynamicSystem, V_List] := MakeDynamicSystem[
	If[MatrixMQ[sys], projectMatrix[MatrixM[sys], V], Null],
	If[MatrixEQ[sys] || FirstOrderSystemQ[sys], projectMatrix[MatrixE[sys], V], Null],
	projectMatrix[MatrixK[sys], V],
	Conjugate[V].MatrixB[sys],
	MatrixC[sys].Transpose[V],
	OutputNames[sys],
	Imtek`Post4MOR`Private`originalMatrixD[sys]
]

Options[ModelReduction] = {ExpansionPoint->0};

ModelReduction[sys_DynamicSystem, dim_Integer, ops___] := Module[
{expnations, start, factor, fun},
	If[Dimensions[MatrixB[sys]][[2]] > 1, Throw["Only a single input is
supported"]];
	expansion = ExpansionPoint /. {ops} /. Options[ModelReduction];
	factor = LinearSolve[(MatrixK[sys] + expansion*MatrixE[sys]), FilterOptions[LinearSolve, ops]];
	start = factor[Flatten[Normal[MatrixB[sys]]]];
	fun = factor[MatrixE[sys].#]&;
	ProjectSystem[sys, ArnoldiProcess[fun, start, dim]]
] /; FirstOrderSystemQ[sys]

ModelReduction[sys_DynamicSystem, dim_Integer, ops___] := Module[
{expnations, start, factor, fun},
	If[Dimensions[MatrixB[sys]][[2]] > 1, Throw["Only a single input is
supported"]];
	expansion = ExpansionPoint /. {ops} /. Options[ModelReduction];
	factor = LinearSolve[(MatrixK[sys] + expansion*MatrixM[sys]), FilterOptions[LinearSolve, ops]];
	start = factor[Flatten[Normal[MatrixB[sys]]]];
	fun = factor[MatrixM[sys].#]&;
	ProjectSystem[sys, ArnoldiProcess[fun, start, dim]]
] /; SecondOrderSystemQ[sys]

Orthogonalize[Vini_List, V_List, tol_:1.*^-12] := Module[{matrixV = Vini, vectorW, norm},
	For[j = 1, j <= Length[V], ++j,
		vectorW = V[[j]];
		For[i = 1, i <= Length[matrixV], ++i,
			vectorW = vectorW - (Conjugate[matrixV[[i]]].vectorW)*matrixV[[i]]
		];
		If[(norm = Norm[vectorW]) < tol,
			Print["Vector is deflated: ", j, " ", norm],
			AppendTo[matrixV, vectorW/norm]
		]
	];
	matrixV
] /; NumericQ[tol]

MakeParametricSystem[matrices_List, matB_, matC_, out_List, fun_List] := ParametricSystem[matrices, matB, matC, out, fun]

Matrices[psys_ParametricSystem] := psys[[1]]

Matrix[psys_ParametricSystem, i_Integer] := psys[[1, i]]

MatrixB[psys_ParametricSystem] := psys[[2]]

MatrixC[psys_ParametricSystem] := psys[[3]]

OutputNames[psys_ParametricSystem] := psys[[4]]

TransformFunctions[psys_ParametricSystem] := psys[[5]]

MatrixMFunction[psys_ParametricSystem] := psys[[5, 1]]

MatrixEFunction[psys_ParametricSystem] := psys[[5, 2]]

MatrixKFunction[psys_ParametricSystem] := psys[[5, 3]]

MakeDynamicSystem[psys_ParametricSystem, par__] := MakeDynamicSystem[psys, {par}]/; VectorQ[{par}, NumericQ]

MakeDynamicSystem[psys_ParametricSystem, par_List] := MakeDynamicSystem[
	MatrixMFunction[psys][Matrices[psys], par],
	MatrixEFunction[psys][Matrices[psys], par],
	MatrixKFunction[psys][Matrices[psys], par],
	MatrixB[psys],
	MatrixC[psys],
	OutputNames[psys]
]

ProjectSystem[psys_ParametricSystem, V_List] := MakeParametricSystem[
	Map[projectMatrix[#, V]& , Matrices[psys]],
	Conjugate[V].MatrixB[psys],
	MatrixC[psys].Transpose[V],
	OutputNames[psys],
	TransformFunctions[psys]
]

End[];

EndPackage[];