(* :Title: Realization *)

(* :Context: Realization` *)

(* :Author: Maris Tonso *)

(* :Summary: This package exports functions to find reduced form of
    discrete-time nonlinear input/output system, state-space realization of
    the input/output system and classical state representation of the
    generalized state-space model. To test if it is possible to perform
    transformations above and to carry out those there are additional
    functions to compute certain sequences of subspaces. *)

(* :Keywords: Nonlinear systems, discrete-time systems, realization, 
    generalized state transformation, input-output models, equivalence,
    linear algebraic method *)

(* :Package Version: 1.0 *)

(* :Mathematica Version: 4.0 *)

(* :Copyright: Copyright 2000 *)

(* :History: *)

(* :Sources:
   1. U.Kotta, P.Liu, A.Zinober. "State-space realization of input-
      output nonlinear difference equations. Proc. European Control
      Conference, Brussels, 1997.

   2. U.Kotta, P.Liu, A.Zinober. "Transfer equivalence of nonlinear
      higher order difference equations". Proc. IFAC Symposium on Nonlinear
      Control Systems Design (NOLCOS), Enschede, 1998, 71-76.

   3. U.Kotta. "Lowering the orderes of input shifts in discrete-time
      generalized state-space systems". Proc. IFAC Conf. on System Structure
      and Control, Nantes, 1998, 759-764.

   4. U. Kotta. Lowering the orders of input shifts in multi-input discrete-
      time generalized state-space systems. Proc. the 37th IEEE Conf. on
      Decision and Control, Tampa, Florida,1998. *)
		
(* :Warnings: The fucntion IntegratePfaff does not work in the case of more
	complecated examples and this may cause the other functions resisting.
	You should realize this is non-commercial, rather experimental software, 
	therefore do not be disappointed if everything does not work as you 
	expect.*)


BeginPackage["Realization`"]



(* -------------------------------------------------------------------- *)
(* ========================== USAGE MESSAGES ========================== *)
(* -------------------------------------------------------------------- *)


(* === Functions, necessary to compute the subspaces {Ik} and {Hk}. === *)

Span::usage =
"Span[{vec1, vec2, ...}, {x1, x2, ... }] represents the difference
vector subspace spanned over the vectors vec1, vec2, ... , where {x1, x2, ...} 
is a list of variables. Example:\n
Span[{{1, x2, 0}, {0, x3, 2}},{x1, x2, x3}]\n"

\[EmptySet]::usage =
"\[EmptySet] is used to represent empty difference subspace 
Span[{},{x1, x2,...}] in TraditionalForm and DifferentialForm."

DifferentialForm::usage =
"DifferentialForm prints vector subspace using differentials. 
DifferentialForm[ Span[{{a1, a2,..}, {b1, b2, ...}, ...}, {x1, x2, ... }] 
---> Span[a1 \[DifferentialD]x1 + a2 \[DifferentialD]x2 + ... , b1
\[DifferentialD]x1 + b2 \[DifferentialD]x2 + ... ].\n
DifferentialForm acts as 'wrapper', whitch affects printing, but not 
evaluation."

TableForm::usage =
"TableForm[list] prints with the elements of list arranged in an array of 
rectangular cells.\n
TableForm[Span[{{a1, a2,..}, {b1, b2, ...}, ...}, {x1, x2, ... }]] 
prints Span-object in form\n
\[DifferentialD]x1  \[DifferentialD]x2  ...\n
 a1    a2  ...\n
 b1    b2  ...\n
 .............\n
TableForm acts as 'wrapper', whitch affects printing, but not evaluation."

Basis::usage =
"Basis[ extEq, {{x1, x2,...}, u}, t ] gives the basis of a vector difference 
subspace, connected with extended state-space system extEq.\n
{x1, x2, .. } - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u."

ForwardShift::usage =
"ForwardShift[ sp, extEq, {{x1, x2,...}, u}, t ] computes the forwardshifted 
set of one-forms from the input set of one-forms.\n
sp - difference vector subspace which represents the input set of one-
forms and must have the form Span[vecs, {x1, x2,..}].\n
extEq - extended state-space system\n
{x1, x2, .. } - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use list of input variables {u1, u2, ...} \n
instead of a single variable u."

BackwardShift::usage =
"BackwardShift[ sp, extEq, {{x1, x2,...}, u}, t ] computes the backwardshifted 
set of one-forms from the input set of one-forms.\n
sp - difference vector subspace which represents the input set of one-
forms and must have the form Span[vecs, {x1, x2,..}]\n
extEq - extended state-space system\n
{x1, x2, .. } - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
BackwardShift[ sp, extEq, {{x1, x2,...}, u}, t , psi]  or \n
BackwardShift[ sp, extEq, {{x1, x2,...}, u}, t , {psi1,psi2,...}] determines 
an operator psi used to compute backwardshift of the set of one-forms.\n 
Alternative operators give the same backward-shifted subspace as a result, but 
the bases may be different. Function BackwardShiftOperator finds all 
alteranative operators psi1, psi2, ... ."

BackwardShiftOperator ::usage =
"BackwardShiftOperator [ extEq, {{x1, x2,...}, u}, t ] finds all alternative 
operators to compute backward shift of a subspace. Alternative operators give 
the same backward-shifted subspace as a result, but the bases may be different.\n
{x1, x2, .. } - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u."

SimplifyBasis::usage =
"SimplifyBasis[sp] simplifies the basis of the vector space sp. The 
space sp must have the form Span[{vec1,vec2,..}, {x1,x2,...}]"

IntersectionSpace::usage =
"IntersectionSpace[sp1, sp2] finds intersection of two vector spaces sp1 
and sp2. Vector spaces must have the form Span[{vec1,vec2,..}, {x1,x2,...}]."

SequenceI::usage =
"SequenceI[ extEq, {{x1, x2, ...}, u }, t] computes the sequence of subspaces
{I1,...,I\[Infinity]} for the state-space equations extEq.\n
{x1, x2, .. } - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
SequenceI[ extEq, {{x1, x2, ...}, u }, t, n ] computes the subspaces {I1,...,In}"

SequenceH::usage =
"SequenceH[ extEq, {{x1, x2, ...}, u }, t] computes the sequence of subspaces 
{H1,...,H\[Infinity]} for the state-space equations extEq.\n
{x1, x2, .. } - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
SequenceH[ extEq, {{x1, x2, ...}, u }, t, n ] computes the subspaces {H1,...,Hn}"

SubmersionQ::usage =
"Submersion[ extEq, {{x1,x2,..}, u}, t] tests if the state-space equations
extEq are submersive or not.\n
{x1,x2,...} - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n"

SameQ::usage =
"SameQ[sp1, sp2] tests if two basis, sp1 and sp2, span one and the same
space or not."

VerifySequence::usage =
"VerifySequence[ extEq, {{x1,x2,...},v}, t] verifies, if the computations of 
the sequences {Ik} and {Hk} for the extended state-space system extEq are correct.
Correct calculations returns list {{True,Null},{True,True},{True,True},...}.\n
{x1,x2,...} - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
Options[VerifySequence] = {ShowInfo}\n
ShowInfo has default value True, which prints information about verified 
conditions.\n
VerifySequence[ extEq, {{x1,x2,...},v}, t, n] tests conditions for first n
subspaces."




(* ============== Exterior differentials and integration =============== *)

Wedge::usage =
"Wedge[df1, df2,...] computes wedge product (exterior product) of the 
differential forms df1, df2,... "

De::usage =
"De[df] computes exterior differential of the differential form df.\n
De[df, Coordinates -> {x1,x2,...}] computes exterior differntial of the 
differential form df, treating symbols x1, x2 as variables.\n

De[x] represents exterior differential dx. In TraditionalForm De[x] is
printed in a form \[DifferentialD]x.\n
De[x,y,z] represents wedge product of differentials dx\[Wedge]dy\[Wedge]dz."

Coordinates::usage =
"Coordinates is an option for De."

ToDifferential::usage =
"ToDifferntial converts Span-object to list of differential one-forms. 
ToDifferential[ Span[{{a1, a2,..}, {b1, b2, ...}, ...}, {x1, x2, ... }]] 
---> { a1De[x1] + a2De[x2] + ..., b1De[x1] + b2De[x2] + ..., ... ]"

Integrability::usage =
"Integrability[sp] tests if the difference subspace of one-forms sp is 
completely integrable or not. This is done by checking the Frobenius 
condition."

IntegratePfaff::usage =
"IntegratePfaff[ {df1,df2,..}, {x1,x2,...} ] tries to integrate Pfaffian 
system of one-forms df1,df2,... treating x1,x2,... as variables.\n
IntegratePfaff[Span[{vec1,vec2,...},{x1,x2,...}]] computes one-forms from 
vector difference subspace and then tries to integrate.\n
General integration of one-forms is a very complicated problem, which has 
no general solution neither a algorithm. The program implemented here works 
only in the most simple cases. If IntegratePfaff fails, try to transform the 
one-forms to exact one-forms and then use IntegratePfaff."




(* ================= Transforming input/output equations ================ *)


IOToExtendedState::usage =
"IOToExtendedState[ ioEq, {y, u}, t] gives the extended state space system for the
input-output equation ioEq.\n
y - output variable\n
u - input variable\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
Options[IOEToESE] = {NewVariables}\n
NewVariables -> {{\!\(z\_1\),\!\(z\_2\),...},{\!\(v\_1\),\!\(v\_2\),...}}
specifies the symbols, used in the extended system for state and input variables.
Deafult value is Automatic, which produces new state variables z1,z2,... and
input variables v1,v2,..."

NewVariables::usage =
"NewVariables is an option for IOEToESE, GSEToESE, IOEToGSE and Realization."

Irreducibility::usage =
"Irreducibility[ ioEq, {y, u}, t] determines whether the nonlinear higher 
order input-output difference equation ioEq is irreducible or not.\n
y - output variable\n
u - input varible\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
Method->SubSpaceH (default value) or Method->SubSpaceI uses sequences Hk or
Ik respectively for testing."

SubSpaceH::usage =
"SubSpaceH is a value for option Method for Irreducibility, Realizability, Realization, 
ClassicTransformability."

SubSpaceI::usage =
"SubSpaceI is a value for option Method for Irreducibility, Realizability, 
ClassicTransformability."

ReducedDifferentialForm::usage =
"ReducedDifferentialForm[ioEq, {y, u}, t ] finds a reduced differntial form 
for reducible input-output difference equation ioEq.\n
y - output variable\n
u - input varible\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u."

Reduction::usage =
"Reduction[ ioEq, {y, u}, t ] finds the reduced i/o difference equation 
from given input-output equation ioEq.\n
y - output variable\n
u - input varible\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u."

Equivalence::usage =
"Equivalence[ ioEq1, ioEq2, {y, u}, t ] determines whether two nonlinear 
higher order input-output difference equations ioEq1 and ioEq2 are 
transfer equivalent or not.\n
y - output variable\n
u - input varible\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u."

Realizability::usage =
"Realizability[ ioEq, {y, u}, t] determines whether the nonlinear higher 
order input-output difference equation ioEq can be realized in the classical 
state-space form.\n
y - output variable\n
u - input varible\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
Method->SubSpaceH (default value) or Method->SubSpaceI uses sequences Hk or
Ik respectively for testing."

Realization::usage =
"Realization[ ioEq, {y, u}, t ] finds the state variables that allow to 
bring the input-output difference equation ioEq into the classical 
state-space form.\n
y - output variable\n
u - input varible\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
Options[Realization] = { Method, ShowInfo, NewVariables}.\n
Method -> {expr1,expr2,...} takes state variables equal to x1=expr1, 
x2=expr2,...\n
Method -> Automatic uses simple models whenever it is possible and 
subspace \!\(H\_s+2\) otherwise.\n
Method -> SubSpaceH  uses subspace \!\(H\_s+2\) for realization.\n
NewVariables -> {{\!\(x\_1\),\!\(x\_2\),...}} determines symbols for 
new state variables. Default value Automatic causes symbols x1,x2,... .\n
ShowInfo -> (True or False) If the value of the option is True, then the 
program informs which method was used for realization. Default value for
ShowInfo is False."

NewVariables::usage =
"NewVariables is an option for IOToExtendedState, IOToGeneralizedState, 
GeneralizedStateToExtendedState, Realization."

ShowInfo::usage =
"ShowInfo is an option for Realization and VerifySequence."




(* ============== Transforming generalized state equations ============== *)

IOToGeneralizedState::usage =
"IOToGeneralizedState[ ioEq, {y, u}, t] finds generalized state equations 
from nonlinear higer order input-output difference equation ioEq.\n
y - output variable\n
u - input varible\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
Options[IOToGeneralizedState] = NewVariables\n
NewVariables -> {{\!\(x\_1\),\!\(x\_2\),...}} determines the symbols for 
new state variables. Default value Automatic produces symbols x1,x2,..."

GeneralizedStateToExtendedState::usage =
"GeneralizedStateToExtendedState[ stEq, {{x1,x2,...}, u}, t] gives the 
extended state space system for the generalized state space model stEq.\n
{x1,x2,...} - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use list of input variables {u1, u2, ...} 
instead of a single variable u.\n
Options[GeneralizedStateToExtendedState] = {NewVariables}\n
NewVariables -> {{\!\(z\_1\),\!\(z\_2\),...},{\!\(v\_1\),\!\(v\_2\),...}}
specifies the symbols, used in the extended system for state and input variables.
Deafult value is Automatic, which produces new state variables z1,z2,... and
input variables v1,v2,..."

Lower::usage =
"Lower[ stEq, {{x1,x2,...}, u}, t] tries to lower the order of the input 
time-shift as much as possible.\n
Options[Lower] = {LoweringOrder, Method}\n
LoweringOrder->n tries to lower the order of the input time-shfit exactly by n.\n
Method -> {expr1,expr2,...} immediately uses the given expressions for 
generalized state transformation."

LoweringOrder::usage = "LoweringOrder is an option for Lower."

ClassicTransformability::usage =
"ClassicTransformability[ stEq, {{x1,x2,...}, u}, t] determines whether the 
generalized state equations stEq can be transformed into the classical 
state-space form.\n
{x1,x2,...} - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
Method->SubSpaceH (default value) or Method->SubSpaceI uses sequences Hk or
Ik, respectively, for testing."

ClassicTransform::usage =
"ClassicTransform[ stEq, {{x1,x2,...}, u}, t] finds the generalized state 
transformation that allows to bring the generalized state equations stEq into
the classical state space form.\n
{x1,x2,...} - state variables\n
u - input variable\n
t - time\n
In case of multi-input system use the list of input variables {u1, u2, ...} 
instead of a single variable u.\n
Options[ClassicTransform] = {Method}\n
Method -> {expr1,expr2,...} immediately uses the given expressions for 
generalized state transformation."








(* -------------------------------------------------------------------- *)
(* ============================= PROGRAMS ============================= *)
(* -------------------------------------------------------------------- *)


Begin["`Private`"]

Off[ General::spell1 ]
epsilon = 10.^-15




(* -------------------------------------------------------------------- *)
(* === Functions, necessary to compute the subspaces {Ik} and {Hk} ==== *)
(* -------------------------------------------------------------------- *)


DifferentialForm[ Span[ sp_?MatrixQ, coords_?VectorQ ]] := Format[
	DisplayForm[ 
		ToBoxes[ Span @@ Dot[ sp, De /@ coords ] ] /.
			{ a___, RowBox[{"De", x___ }], " ", b__ } ->
		{ a , b, " ", RowBox[{ "De", x }] } /.
			{ "De", "[", x_, "]" } -> { "\[DifferentialD]", x }
	] (* End DisplayForm *)
] (* End Format *)

DifferentialForm[ sps_List ] := DifferentialForm /@ sps
DifferentialForm[ Span[ { }, _] ] := Format[ \[EmptySet] ]
DifferentialForm[ ___ ] := err /; Message[ DifferentialForm::usage ]



Unprotect[ TableForm ];
TableForm[ Span[ sp_?MatrixQ, coords_?VectorQ ]] := 
	TableForm[ sp, 
                TableHeadings -> { None,
                        "\[DifferentialD]"<>ToString[#]& /@ coords
                },
		TableAlignments -> Center
	] (* End TableForm *)
TableForm[ Span[ { }, coords_] ] := Format[ \[EmptySet] ]



Format[ Span[ sp_?MatrixQ, coords_?VectorQ ], TraditionalForm ] := 
DisplayForm[
	ToBoxes[ sp.De /@ coords, TraditionalForm ] /. 
	FormBox[ InterpretationBox[
                RowBox[{"\[DifferentialD]", rbox_}], __], TraditionalForm] -> 
		RowBox[{"\[DifferentialD]", rbox}] /.
        { RowBox[{"\[DifferentialD]", a___ }], " ", rboxs__ } ->
                { rboxs, " ", RowBox[{ "\[DifferentialD]", a }] } /.
	({#, "[", i_ ,"]" } -> {#, "(", i, ")" }& /@ 
		( ToBoxes[ Head@# ]& /@ coords )) /.
	FormBox[ rbox_, TraditionalForm ] -> 
		FormBox[ RowBox[{"span", rbox} ] , TraditionalForm ]
] (* End DisplayForm *)
Format[ Span[ { }, _], TraditionalForm ] := \[EmptySet]



Basis[ extEq_, {Z_, v_?(Head[#] =!= List&)}, t_] :=
	Basis[ extEq, {Z, {v}}, t]

Basis[ extEq:{(_[t_+1]==_)..}, {Z:{__},V:{__}}, t_] :=
	Span[ Take[ 
		IdentityMatrix[ Length@Z + Length@V ], Length@Z ], 
	#[0]& /@ Join[Z,V] ]

Basis[ ___ ] := err /; Message[ Basis::usage ]



ForwardShift[ Span[ { }, coords_], extEq_, {Z_, V_}, t_] := 
	Span[ { }, coords ]

ForwardShift[ sp_, extEq_, {Z_, v_?(Head[#] =!= List&)}, t_] :=
	ForwardShift[ sp, extEq, {Z,{v}}, t]

ForwardShift[ Span[ sp_, coords_], extEq:{(_[t_+1]==_)..}, 
	{Z:{__}, V:{__}}, t_] := 
With[ { rules = extEq /. t+_. -> 0 /. Equal -> Rule},
Span[
	Dot[ Take[#, Length@Z ]& /@ sp /. 
		Join[ rules, 
			#[j_] -> #[j+1]& /@ V ,
			#[j_?Negative] -> #[j+1]& /@ Z
		],
		Outer[ D, #[0]&/@Z /. rules, coords]
	], 
coords ]]

ForwardShift[ ___ ] := err /; Message[ ForwardShift::usage ]



BackwardShiftOperator[ extEq_, {Z_, v_?(Head[#] =!= List&)}, t_] := 
		BackwardShiftOperator[ extEq, {Z, {v}}, t ]

BackwardShiftOperator[ extEq:{(_[t_+1]==_)..},
	{Z:{__}, V:{__}}, t_] := 
Module[{ dz0, dz1, int, compsp, w0, 
	z0 = #[0]&/@Z }, 
	dz0 = Span[ IdentityMatrix[ Length@Z ], z0];
	dz1 = ForwardShift[ dz0, extEq, {Z, V}, t];
	int = IntersectionSpace[ dz0, dz1]; 
	compsp = NullSpace[ First@int ];
	If[ compsp === {},
		Solve[ extEq /. t->0, z0, #[0]& /@ V ], 
	(* Else *)
		w0 = z0[[#]]& /@ Distribute[ 
			Flatten[ Position[ #, _?(Chop[#] =!= 0&), 
				{1}, Heads -> False]]& /@ compsp,
		List ];
		Flatten[ Solve[ 
			extEq /. t->0, Complement[z0, #], #[0]& /@ V
		]& /@ w0, 1]
	] (* End If *)
] (* End Module *)

BackwardShiftOperator[ ___ ] :=
    err /; Message[ BackwardShiftOperator::usage ]




BackwardShift[ Span[ { }, coords_], __] := Span[ { }, coords ]

BackwardShift[ sp_, extEq_, {Z_, v_?(Head[#] =!= List&)}, t_,
psi___] := BackwardShift[ sp, extEq, {Z,{v}}, t, psi ]

BackwardShift[ sp_, extEq_, {Z_, V_}, t_, psi_?MatrixQ] := 
	BackwardShift[ sp, extEq, {Z, V}, t, #]& /@ psi

BackwardShift[ sp_Span, extEq_, {Z_, V_}, t_ ] := 
Module[{ psi, bwssp},
	psi = BackwardShiftOperator[ extEq, {Z, V}, t];
	bwssp = BackwardShift[ sp, extEq, {Z, V}, t, #]& /@ psi;
	First @ Select[ bwssp, 
		LeafCount[#] == Min[LeafCount /@ bwssp]&, 1 ]
] (* End Module*)

BackwardShift[ Span[ sp_, coords_], extEq:{(_[t_+1]==_)..}, 
	{Z:{__}, V:{__}}, t_, psi_?VectorQ ] := 
Module[ { z0 = #[0]&/@Z, rules },
	rules = MapThread[ Rule, 
		{z0, z0 /. psi /. (#[j_] -> #[j-1]& /@ Z) } 
	];
	Span[
		Dot[ Take[#, Length@Z ]& /@ sp /. 
			Join[ rules, 
				#[j_] -> #[j-1]& /@ V
			], (* End Join *)
			Outer[ D, z0 /. rules, coords ]
		], (* End Dot *) 
	coords ] (* End Span *)
] (* End Module *)

BackwardShift[ ___ ] := err /; Message[ BackwardShift::usage ]



SimplifyBasis[ Span[{ }, coords_]] := Span[ { }, coords ]

SimplifyBasis[ Span[ sp_?MatrixQ, coords_]] := 
Span[
Module[ {fact, vec1,
	sp1 = sp /. c_Real?(Round[#] == #& ) -> Round[c] },
	sp1 = If[ LeafCount[#1] < LeafCount[#2], #1, #2]& 
		[   Chop[ Simplify[sp1], epsilon ],
			Chop[ Simplify[ RowReduce[sp1] ], epsilon ]
		];
	Function[ vec,
		vec1 = PowerExpand[ Expand[
			PolynomialLCM@@
				Denominator[ Together@vec ] * vec
		] ];
	 	ExpandAll[ Simplify[ 
	 		If[ fact = PolynomialGCD@@vec1; fact=!=0, 
	 			vec1/fact, vec1
	 		] (* End If *)
	 	] ]
	 ] /@ sp1
] (* End Module *) /. c_Real?(Round[#] == #& ) -> Round[c],
coords ] (* End Span *)

SimplifyBasis[ ___ ] := err /; Message[ SimplifyBasis::usage ]


IntersectionSpace[ Span[sp1_, coords_], Span[sp2_, coords_] ] := 
Span[ 
With[ {
	compsp = Join[ NullSpace@ sp1, NullSpace@ sp2 ]/.
	c_Real?(Round[#] == #&) -> Round[c]},
	If[ compsp === {},
		IdentityMatrix[ Length@ coords ],
		NullSpace[ compsp ] /. 
			c_Real?(Round[#] == #&) -> Round[c]
	] (* End If *)
], (* End With *)
coords ]

IntersectionSpace[ ___ ] := err /; Message[ IntersectionSpace::usage ]



SequenceI[ extEq:{(_[t_+1]==_)..}, {Z:{__},V_}, t_, n_:10 ] := 
FixedPointList[
    SimplifyBasis[
        IntersectionSpace[#, ForwardShift[#, extEq, {Z, V}, t]]
    ]&,
    Basis[ extEq, {Z, V}, t],
    n - 1, 
    SameTest -> (#1==#2 || #1[[1]]==-#2[[1]] || #2[[1]]=={} & )
] (* End FixedPointList *)

SequenceI[ ___ ] := err /; Message[ SequenceI::usage ]



SequenceH[ extEq:{(_[t_+1]==_)..}, {Z:{__},V_}, t_, n_:10 ] := 
Module[ { H1, H2, intsp, psi },
	H1 = Basis[ extEq, {Z, V}, t];
	intsp = IntersectionSpace[ H1, 
		ForwardShift[ H1, extEq, {Z,V}, t] ];
	psi = BackwardShiftOperator[ extEq, {Z,V}, t];
	H2 = BackwardShift[ intsp, extEq, {Z,V}, t, psi];
	psi = psi[[ First[ First[ 
		Position[ H2, 
			_?(LeafCount[#] == Min[ LeafCount /@ H2] & ),
		 1 ]  
	]]  ]];
	FixedPointList[
		SimplifyBasis[ 
			BackwardShift[
				IntersectionSpace[ #, 
					SimplifyBasis[
						ForwardShift[ #, extEq, {Z,V}, t ]
					] (* End SimplifyBasis *)
				], (* End IntersectionSpace *)
			extEq, {Z,V}, t, psi]
		]&,
		H1, n - 1, 
		SameTest -> (#1==#2 || #1[[1]]==-#2[[1]] || #2[[1]]=={} & )
	] (* End FixedPointList *)
] (* End Module *)

SequenceH[ ___ ] := err /; Message[ SequenceH::usage ]



SubmersionQ[ extEq_, {Z_, v_/;Head[v] =!= List}, t_] :=
    SubmersionQ[ extEq, {Z, {v}}, t]

SubmersionQ[ extEq:{(_[t_+1]==_)..}, {Z:{__},V:{__}}, t_] :=
	With[{
        f = extEq /. Equal[_, fi_] -> fi,
		vars = #[t]& /@ Join[Z,V] },
		Count[ 
			RowReduce@ Transpose@ Outer[ D, f, vars ],
			{___?(#==0&), 1, ___?(#==0&)}
		] == Length@f
	] (* End With *)

SubmersionQ[ ___ ] := err /; Message[ SubmersionQ::usage ]



Unprotect[ SameQ ];
SameQ[ Span[ sp1_, coords_], Span[ sp2_, coords_] ] :=
If[ sp1 === sp2, 
	True,
	Apply[ SameQ, 
		Simplify[ RowReduce[#] //.
			{vecs___, {___?(#==0&)} } -> {vecs} 
		]& /@ {sp1, sp2}
	] (* End Apply *)
] (* End If *) 



VerifySequence[ extEq:{(_[t_ + 1] == _)..}, {Z:{__}, V_}, t_, 
	n_:10, opts___] := 
Module[ {Ik, Hk, Istar, shi, Hiplus, res }, 
	Ik = SequenceI[extEq, {Z, V}, t, n]; 
	Hk = SequenceH[extEq, {Z, V}, t, n];
    shi = ShowInfo /. {opts} /. Options[ VerifySequence ];
	Table[
		Istar = Nest[ 
			SimplifyBasis[ ForwardShift[#, extEq, {Z,V}, t]]& , 
		Hk[[i]], i - 1 ];
		Hiplus = SimplifyBasis[
			ForwardShift[ Hk[[i]], extEq, {Z, V}, t]
		];		
		res = { SameQ[ Istar, Ik[[i]] ],
			If[ i > 1,
				SameQ[ Hiplus,
					IntersectionSpace[ Hiplus, Hk[[i-1]] ]
				] (* End SameQ*)
			] (* End If *)	
		};
		If[ shi,
			Print["I", i, " = ", Ik[[i]] ]; 
			Print["H", i, " = ", Hk[[i]] ];
			Print[ "I", i, "\[FivePointedStar] = \[CapitalDelta]", 
            	Superscript["", i - 1], "(H", i, ") = ", Istar];
            Print[ "I", i, If[ res[[1]], " = ",
                "!= "], "I", i, "\[FivePointedStar]"];
            If[ i > 1,
	        	Print["H", Superscript[ i, "+"], " = ", Hiplus ];
 	        	Print["H", Superscript[ i, "+"], 
  	        		If[ res[[2]], " \[Element] ", " \[NotElement] " ], 
            	"H", i-1];
            ];
            Print[]
          
		];
		res,
	{ i, Length@Hk } ] (* EndTable *)
] (* End Module*)

VerifySequence[ ___ ] := err /; Message[ VerifySequence::usage ]

Options[ VerifySequence ] = { ShowInfo -> True }






(* --------------------------------------------------------------------- *)
(* ============== Exterior differentials and integration =============== *)
(* --------------------------------------------------------------------- *)

De[x__?(Head[#]=!=Rule&)] := Signature[{x}] De @@ Sort[{x}] /; 
	!OrderedQ[{x}]
De[x__] := 0 /; Signature[{x}] == 0

De[ df:A_. De[x__], opts___Rule]:=
Module[{ coords = Coordinates /.
		{opts} /. Options[ De ]},
			If[ coords === Automatic,
				coords = Union[ Cases[ df, _Symbol, -1]]
			];
       Wedge[ Plus @@ ( D[A, #] De[#]& /@ coords), De[x] ] 
]

De[A_. df_Plus, opts___Rule] :=
	     De[#, opts]& /@ Expand[A df] /; !FreeQ[A df, De ]

(*De[ ___ ] := err /; Message[ De::usage ]*)


De /: MakeBoxes[ De[x_, y___], TraditionalForm ] := 
RowBox[{ "\[DifferentialD]", MakeBoxes[x, TraditionalForm ],
	If[ Length[{y}] === 0, "", 
		seq[ "\[Wedge]", MakeBoxes[ De[y], TraditionalForm ]]
	] (* End If *)
}] /. seq -> Sequence


Options[ De ] = {Coordinates -> Automatic}



Wedge[df_] := df

Wedge[A_. De[x__], B_. De[y__], df___] := Wedge[ A B De[x,y], df]

Wedge[df1___, f_, df2___] := f Wedge[ df1, df2 ] /; 
	FreeQ[f, De[__]]

Wedge[df1___, A_. df2_Plus, df3___] := 
	Wedge[ df1, Expand[A #], df3]& /@ df2



ToDifferential[ Span[sp_?MatrixQ, coords_?VectorQ]] := Dot[ sp, De /@ coords ]
ToDifferential[ Span[{}, _?VectorQ ]] := {}
ToDifferential[ sps_List] := ToDifferential /@ sps
ToDifferential[ ___ ] := err /; Message[ ToDifferential::usage ]



Integrability[ Span[{ }, _] ] := {}

Integrability[ Span[ sp_, coords_ ]] := 
Module[{ pforms = Dot[sp, De/@coords ], coords1 },
coords1 = Union@ Cases[ pforms, Alternatives @@ 
    (coords /. x_[0] -> x[_]), -1];
If[ Complement[ coords1, coords] =!= {},
	False,
(* Else *)
	Apply[ And, 
		Chop[
			Simplify[
				Wedge[
					De[ #, Coordinates -> coords1 ],
				Wedge @@ pforms ] (* End Wedge *)
			], (* End Simplify *)
		epsilon ] (* End Chop *)
	=== 0& /@ pforms ] (* End Apply *)
] (* End If *)
] (* End Module *)


Integrability[ ___ ] := err /; Message[ Integrability::usage ]



IntegratingFactor[ A1_. De[x1_] + A2_. De[x2_] ] :=
Module[ { my, x, y },
	First[ my[x,y] /. DSolve[
		D[ (A1 /. {x1->x, x2->y}) my[x,y], y] == 
		D[ (A2 /. {x1->x, x2->y}) my[x,y], x],
	my[x,y], {x,y} ] ] /.
	{x->x1, y->x2}
] (* End Module *)

IntegratingFactor[ __ ] := $Failed




IntegratePfaff::unable = "IntegratePfaff is unable to transform 
the set of one-forms to an exact one-forms.";

IntegratePfaff::nonint = "The subspace is not completely integrable. "

IntegratePfaff[ Span[ { }, _] ] := { }

IntegratePfaff[ Span[ sp_, coords_] ] := 
	IntegratePfaff[ ToDifferential@Span[ sp, coords ], coords ]
	
IntegratePfaff[ dforms:{__}, coords0_] :=
Module[{ 
	dforms1 = Expand /@ dforms,
	coords = Union@ Cases[ dforms, 
		Alternatives @@ (coords0 /. x_[0] -> x[_]), -1 ],
	dfnew, ddf, wddf, combin, dfothers, df1, df1list, df2, k,
	len, M, A, x, B, y, totdif, omega, beeta, mu, eqv, sol, F, s, result },
result = Function[ df,
dfnew = Which[ 
	(ddf = Chop[ De[ df, Coordinates -> coords ], epsilon ]) === 0,
    	df,
    Chop[ Wedge[ ddf, Wedge@@dforms1 ], epsilon ] =!= 0,
    	nonintegrable,
    (wddf = Chop[ Wedge[ ddf, df ], epsilon ]) === 0,
    	mu = IntegratingFactor[ df ] /. C[1][_] -> 1; 
		If[ mu =!= $Failed, df mu, unable ],
    True,	
		dfothers = Complement[ dforms1, {df} ];
		combin = Distribute[ {{}, {#}}& /@ Range[ Length@dfothers ], 
			List, List, List, Union 
		] // Sort // Rest;
		dfothers = First@Select[ dfothers[[#]]& /@ combin,
			Chop[ Wedge[ wddf, Wedge@@# ], epsilon ] === 0&, 
			1 ]; (* End Select *)
		df1 = df; (* siin on Head[df]===Plus kindlasti *)
		df2 = 0;
		k = 0;
		While[ k++; (* Print[k];
			Print["df1 = ", df1];
			Print["df2 = ", df2];*)
			df1list = toList[df1];
			len = Function[ Adx,
				Length@Select[ coords, !FreeQ[ Adx, #]& ] 
			]/@ df1list;
			M = Max@len;
		M >= 2 && k <= 30,
			{A, x} = df1list[[
				Position[ len, _?(# == M&) ] [[1,1]]
			]] /. A1_ De[x1_] -> {A1, x1};
			(*Print["A De[x] = ", A, "\[DifferentialD]", x ];*)
			y = Union@Complement[
				Select[ coords, !FreeQ[A De[x], #] &], {x} ];
			B = D[ Integrate[A, x], #]& /@ y;
			totdif = A De[x] + B.De/@y;
			Which[
			(* ======== T\[ADoubleDot]isdiferentsiaal on olemas *)
			Length@df1list - Length@toList@Chop[df1-totdif] == M,
				(*Print["Taisdiferentsiaal on olemas"];*)
				df1 = Chop[ df1 - totdif ];
				df2 = df2 + totdif,
			(* ======== Diferentsiaali v\[ADoubleDot]ljataandamine *)
			!FreeQ[ dfothers, De[x] ] && (
				omega = addDifForm[ A, x, dfothers ];
				Length@toList@Chop[ df1 - omega ] < Length@df1list ),
				(*Print["Diferentsiaali v\[ADoubleDot]ljataandamine"];*)
				df1 = Chop[ df1 - omega ],
			(* ======== Teiste komponentide juurdeliitmine *)
			omega = Inner[
				beeta = Coefficient[ df1, De[#2] ];
				Which[
				Chop[ beeta - #1 ] === 0,
					0,
				!FreeQ[ dfothers, De[#2] ],
					addDifForm[ #1-beeta, #2, dfothers ],
				True,
					False
				]&, (* End Which *)		
			B, y ]; (* End Inner *)
			FreeQ[ omega, False],
				(*Print[" Teiste komponentide juurdeliitmine"];*)
				df1 += omega;
				df1 = Chop[ df1 - totdif ];
				df2 += totdif,
			(* ======== Osaline integreeruvustegur *)
			Length@y === 1 && (
			omega = Select[ Cases[ df1list, _. De[ First@y ] ],
				Chop@ Wedge[ De[ A De[x] + #, Coordinates -> coords ],
					A De[x] + # ] === 0& ];
			omega =!= {} ),
				(*Print["Osaline integreeruvustegur"];*)
				mu = IntegratingFactor[ First@omega + A De[x] ] /.
					C[1][_] -> 1;
				df1 = Expand[ mu df ];
				df2 = 0
			] (* End Which *)
		]; (* End While *)
		If[ Chop[ De[ df1 + df2, Coordinates -> coords ], epsilon ] === 0, 
			df1 + df2, unable ]
	]; (* End Which *)

If[ dfnew === unable || dfnew === nonintegrable, 
	dfnew,
	eqv = dfnew /. (#->#[s]& /@ coords) /.
		De[zi_[s]] -> zi'[s];
	N@ Expand[ 
		F[s] /. First@
		DSolve[ eqv == F'[s], {F[s]}, s] /. 
		C[1] -> 0 /. zi_[s] -> zi
	] /. c_Real?(Round[#]==#&) -> Round[c]
] (* End If *)

] (* End Function *) /@ dforms1 /. x_Real :> N[x,6];

If[ FreeQ[ result, nonintegrable], 
	result, 
	Message[ IntegratePfaff::nonint ]; { } 
	] 
/; FreeQ[result, unable] ] (* End Module *)

(* Kui integreerimine on vo~imalik, kuid ei oska seda teha.*)
IntegratePfaff[ dforms:{__}, coords0_] := unable /; 
	Message[ IntegratePfaff::unable, df ]

IntegratePfaff[ ___ ] := err /; Message[ IntegratePfaff::usage ]

toList[x_] := Switch[x, _Plus, List@@x, 0, {}, _, {x} ]

addDifForm[ A_, x_, dfothers_] :=
Module[ { omega },				
	omega = Select[ dfothers, !FreeQ[ #, De[x] ]& ];
	omega = Flatten@Cases[ omega, 
		coef_. De[x] -> Expand[ A omega / coef ], {0, 2}];
	First@Select[ omega, 
		LeafCount[#] == Min[ LeafCount /@ omega ]&, 1]
] (* End Module *)







(* ---------------------------------------------------------------------- *)
(* ================= Transforming input/output equations ================ *)
(* ---------------------------------------------------------------------- *)



IOToExtendedState::nvars = "Value of option NewVariables -> `1` should be a list
containing two list of variables: first should have length `2` and the 
second length `3`.";

IOToExtendedState[ ioEq_, {y_, u_/;Head[u]=!=List}, t_ , opts___] :=
	IOToExtendedState[ ioEq, {y, {u}}, t, opts]

IOToExtendedState[ ioEq_Equal, {y_, U_List}, t_ , opts___] :=
Module[ {n, m, s, nvars, X, stEq, extEq },
	n = Max@ Cases[ ioEq, y[t+i_.] -> i, -1 ];
	m = Length@U;
	s = Max[ Cases[ ioEq, #[t+i_.] -> i, -1]& /@ U];
	nvars = NewVariables /. {opts} /. Options[ IOToExtendedState ];
    {Z, V} = If[ Length[ nvars ] == 2 && 
			Length[ nvars[[1]] ] == n + m(s+1) &&
			Length[ nvars[[2]] ] == m, 
		nvars,
	(* Else *)
		If[ nvars =!= Automatic,
			Message[ GeneralizedStateToExtendedState::nvars, 
				nvars, n + m(s+1) , m];
		]; (* End If *)
		{ Array[ ToExpression["z" <> ToString[#] ]&, n + m(s+1) ],
		  Array[ ToExpression["v" <> ToString[#] ]&, m ] }
	]; (* End If *)
	X = Array[ ToExpression[ "x" <> ToString[#]]&, n];
	stEq = IOToGeneralizedState[ ioEq, {y, U}, t ];
	extEq = GeneralizedStateToExtendedState[ stEq[[1]], {X, U}, t , 
		NewVariables -> {Z, V} ];
	{ extEq[[1]], extEq[[2]] /. Reverse /@ stEq[[2]] }
]

IOToExtendedState[ ___ ] := err /; Message[ IOToExtendedState::usage ]

Options[ IOToExtendedState ] = { NewVariables -> Automatic }



Irreducibility[ ioEq_, {y_, u_/;Head[u]=!=List}, t_, opts___] :=
    Irreducibility[ ioEq, {y, {u}}, t, opts ]

Irreducibility[ ioEq_Equal, {y_/;Head[y]=!=List, U_List}, t_, opts___] :=
Module[ {extEq, rules, Z, V},
    {extEq, rules} = IOToExtendedState[ ioEq, {y, U}, t];
	Z = #[[2,0]]& /@ rules;
    V = Array[ ToExpression[ "v" <> ToString[#]]&, Length@U ];
    met = Method /. {opts} /. Options[ Irreducibility ];
	First @ Last @ #[ extEq, {Z, V}, t ] === {} & 
		[ If[ met === SubSpaceI, SequenceI, SequenceH ] ]
] (* End Module *)

Irreducibility[ ___ ] := err /; Message[ Irreducibility::usage ]

Options[ Irreducibility ] = { Method->SubSpaceH }



ReducedDifferentialForm::irreducibility = "The system `1` is
irreducible and does not have a reduced differential form. ";

ReducedDifferentialForm[ ioEq_, {y_, u_/;Head[u]=!=List}, t_] :=
    ReducedDifferentialForm[ ioEq, {y, {u}}, t ]

ReducedDifferentialForm[ ioEq_Equal, {y_/;Head[y]=!=List, U_List}, t_] :=
Module[ { n, extEq, rules, dphi }, 
	n = Max[ Cases[ ioEq, y[t + i_.] -> i, -1]];
    {extEq, rules} = IOToExtendedState[ ioEq, {y, U}, t];
    Z = #[[2,0]]& /@ rules;
    V = Array[ ToExpression[ "v" <> ToString[#]]&, Length@U ]; 
    dphi = Last @ SequenceH[ extEq, {Z, V}, t] /.
    	(#[0] -> #[t]& /@ Z) /. Reverse /@ rules;
    If[ First[dphi] === { },
		Message[ ReducedDifferentialForm::irreducibility, ioEq ];
		{ },
		dphi 
	] (* End If *)
] (* End Module *)

ReducedDifferentialForm[ ___ ] :=
    err /; Message[ ReducedDifferentialForm::usage ]



Reduction[ ioEq_Equal, {y_/;Head[y]=!=List, U_}, t_] := 
Module[{dphi, ioEq1 = ioEq}, 
    Off[ ReducedDifferentialForm::irreducibility ]; 
    While[
        dphi = ReducedDifferentialForm[ ioEq1, {y, U}, t ];
        dphi =!= {},         
        ioEq1 = First[ IntegratePfaff[ dphi ] ] == 0
    ];
    On[ ReducedDifferentialForm::irreducibility ];
ioEq1 ]

Reduction[ ___ ] := err /; Message[ Reduction::usage ]




Equivalence[ ioEq1_Equal, ioEq2_Equal, 
	{y_ /; Head[y] =!= List, u_}, t_] := 
With[{
	reduced1 = (#1 - #2&) @@ Reduction[ ioEq1, {y, u}, t], 
    reduced2 = (#1 - #2&) @@ Reduction[ ioEq2, {y, u}, t] }, 
    reduced1 === reduced2 || reduced1 === -reduced2
] (* End With *)

Equivalence[ ___ ] := err /; Message[ Equivalence::usage ]



Realizability[ ioEq_Equal, { y_, u_ /; Head[u]=!=List}, t_ , opts___] :=
        Realizability[ ioEq,{y, {u}}, t, opts ]

Realizability[ ioEq_Equal, { y_ /; Head[y]=!=List, U_List}, t_, opts___ ] := 
Module[ {n, m, s, met, extEq, Z, V }, 
    n = Max[ Cases[ ioEq, y[t+i_.] -> i, -1 ]]; 
    s = Max[ Cases[ ioEq, #[t+i_.] -> i, -1]& /@ U];
    m = Length @U;
    met = Method /. {opts} /. Options[ Realizability ];
    extEq = First@ IOToExtendedState[ ioEq, {y, U}, t];
    Z = Array[ ToExpression[ "z" <> ToString[#]]&, n+s+1 ];
    V = Array[ ToExpression[ "v" <> ToString[#]]&, m ];
	Integrability[ Last[
		#[ extEq, {Z, V}, t, s+2 ]
	]] & [ If[ met === SubSpaceI, SequenceI, SequenceH ] ]
]

Realizability[ ___ ] := err /; Message[ Realizability::usage ]

Options[ Realizability ] = { Method -> SubSpaceH }



Realization::method = "Realization uses `1`.";

Realization[ ioEq_, {y_, u_ /; Head[u]=!=List}, t_, opts___] :=
	Realization[ ioEq, {y, {u}}, t, opts]
	
Realization[ ioEq_Equal, {y_/;Head[y]=!=List, U:{__}}, 
	t_, opts___Rule] := 
Module[{ met, shi, n, s, m, X, Z, V, extEq, rules, phi, classic, stvars, 
	classicinvrule,
	u, k, f, shiftmin, shiftmax, result, unabletointgrate },

    met = Method /. {opts} /. Options[ Realization ];
    shi = ShowInfo /. {opts} /. Options[ Realization ];
    X = NewVariables /. {opts} /. Options[ Realization ];
    n = Max[ Cases[ ioEq, y[t+i_.] -> i, -1 ]]; 
    s = Max[ Cases[ ioEq, #[t+i_.] -> i, -1]& /@ U];
    m = Length@U;
    X = If[ Dimensions[X] == {1, n},
		First[X],
	(* Else *)
		If[ X =!= Automatic,
			Message[ IOToGeneralizedState::nvars, X, n ];
		]; (* End If *)
		Array[ ToExpression[ "x" <> ToString[#]]&, n]
	]; (* End If *)
    phi = First@ First@ Solve[ ioEq, y[t+n] ];
    

Which[
	(* ********************** simple models type 1 ************************* *)
	And[ met === Automatic, m == 1,
		u = First@U;
		k = n - s - 1;
		f = Array[ 0&, n - k ];
		result = Catch[ Scan[ (
			{shiftmin, shiftmax} = {Min[#], Max[#]} & 
				[Cases[#, u[t+i_.] | y[t+i_.] -> i, {0, -1 }]];
			If[ shiftmax - shiftmin <= k &&
				If[ FreeQ[#, u ], True,
					Equal@@# && #[[1]] == shiftmin & 
						[Cases[ #, u[t+i_.] -> i, {0, -1 }]]
				], (* End If *)
			(* Then *)
				f[[n - shiftmin - k]] += # /. 
					u[_] -> u[t] /. 
					y[t + i_.] :> X[[i + 1 - shiftmin ]][t],
			(* Else *)
				Throw[ $Failed ]
			] (* End If *)
			)& , toList[ Expand //@ Last[phi] ] 
		]]; (* End Scan & Catch *) 
		result =!= $Failed
	], (* End And  ------------------------------------------------------------- *)
	If[ shi, Message[ Realization::method,"simple models method" ]];
	f = Join[ Array[ 0&, k ], f ];
	classic = MapThread[ #1[t+1] == #2[t] + #3&, 
		{ X, RotateLeft@X, f } ];
	classic[[-1, 2]] -= First[X][t],

	
	(* ********************** simple models type 2 ************************** *)

	(* ********************** given state variables ************************* *)
	And[ Head[met] === List, Length[met] == n ],

	If[ shi, Message[ Realization::method, "given variables" ]];
	stvars = met,
	
	(* *************************** SubSpace Hs+2 ****************************** *)
	True,
		{extEq, rules} = IOToExtendedState[ ioEq, {y, U}, t ];
    	Z = Array[ ToExpression[ "z" <> ToString[#]]&, n + m(s+1) ];
		V = Array[ ToExpression[ "v" <> ToString[#]]&, m ];
    	stvars = IntegratePfaff[ 
    		Last@ SequenceH[ extEq, {Z, V}, t, s+2] ];
		If[ FreeQ[ stvars, De ], 
			stvars = stvars /. ((#[0] -> #[t]& /@ Z) /. Reverse /@ rules)
		]
]; (* End Which *)

	(* *************************** Main program ****************************** *)
result = Which[
	ValueQ[classic],
		Append[ classic, First[X][t] == y[t] ],
	stvars === {},
		{},
	!FreeQ[stvars, De],
		unabletointgrate = 1,
	True,
		If[ shi, Message[ Realization::method, "subspace Hs+2 for realization." ]];
		classic = MapThread[ #1[t] == #2&, { X, stvars }];
		classicinvrule = First@ 
			Solve[ classic, Union@Cases[classic, y[t+_.],-1]  ] /.
			c_Real?(Round[#] == #&) -> Round[c];
		classic = Chop[ 
    		Expand //@ (classic /.
    		t -> t+1 /. phi ), epsilon ] (* End Chop *) /.
			classicinvrule;
		classic = Simplify[ Chop[ Expand //@ classic, epsilon] ] /. 
			c_Real?(Round[#] == #&) -> Round[c];
		Append[ classic, First[X][t] == y[t] ]
];  (* End Which *)
result /; Not[ ValueQ[ unabletointgrate ] ]
] (* End Module*)


Realization[ ___ ] := err /; Message[ Realization::usage ]

Options[ Realization ] = {
	Method -> Automatic,
    ShowInfo -> False, 
	NewVariables -> Automatic 
};



IOToGeneralizedState::nvars = "Value of option NewVariables -> `1` should have a form
{{x1, ... , x`2`}}.";

IOToGeneralizedState[ ioEq_, {y_, u_/;Head[u]=!=List}, t_, opts___] := 
	IOToGeneralizedState[ ioEq, {y, {u}}, t, opts ]

IOToGeneralizedState[ ioEq_, {y_, u_}, t_, opts___]:=
	Module[{n, phi, X , f, rules},
		n = Max@ Cases[ioEq, y[t+i_.] -> i,-1] ;
		phi = y[t+n] /. Solve[ ioEq, y[t+n]] // First;
		X = NewVariables /. {opts} /. Options[ IOToGeneralizedState ];
		X = If[ Dimensions[X] == {1, n},
			First[X],
		(* Else *)
			If[ X =!= Automatic,
				Message[ IOToGeneralizedState::nvars, X, n ];
			]; (* End If *)
			Array[ ToExpression[ "x" <> ToString[#]]&, n]
		]; (* End If *)
		rules = Table[ y[t+i-1] -> X[[i]][t], {i,n} ];
		f = Append[ 
			#[t]& /@ Rest[X], phi ];
		{ MapThread[ Equal, {#[t+1]& /@ X, f}] /. rules,
			rules }
	] (* End Module *)

IOToGeneralizedState[ ___ ] := err /; Message[ IOToGeneralizedState::usage ]
	
Options[ IOToGeneralizedState ] = { NewVariables -> Automatic }




(* ---------------------------------------------------------------------- *)
(* ============== Transforming generalized state equations ============== *)
(* ---------------------------------------------------------------------- *)


GeneralizedStateToExtendedState::nvars = "Value of option NewVariables -> `1` 
should be a list containing two list of variables: first should have length `2` 
and the second length `3`.";

GeneralizedStateToExtendedState[ stEq_, {X_, u_/;Head[u]=!=List}, t_, opts___] :=
    GeneralizedStateToExtendedState[ stEq, {X, {u}}, t, opts ]
	
GeneralizedStateToExtendedState[ stEq:{(_[t_+1]==_)..}, {X:{__}, U:{__}}, t_, opts___] := 
	Module[ {n, alpha, m, Z, V, extEq, l, rules, rules1, nvars },
	n = Length@X; 
	alpha = Max[ Cases[ stEq, #[t+i_.]->i, -1]& /@ U];
	m = Length@U;
	nvars = NewVariables /. {opts} /. Options[ GeneralizedStateToExtendedState ];
	{Z, V} = If[ Dimensions[ nvars ] == {2} && 
			Length[ nvars[[1]] ] == n + m(alpha+1) &&
			Length[ nvars[[2]] ] == m,
		nvars,
	(* Else *)
		If[ nvars =!= Automatic,
			Message[ GeneralizedStateToExtendedState::nvars, nvars, n + m(alpha+1) , m];
		]; (* End If *)
		{ Array[ ToExpression[ "z" <> ToString[#]]&, n + m(alpha+1) ],
          Array[ ToExpression[ "v" <> ToString[#]]&, m ]}
	]; (* End If *)
	extEq = stEq;
	rules = MapThread[ #1[t]->#2[t]&, {X, Take[Z,n]} ];
	rules1 = rules /. t -> t+1;
 	Do[ 
		l = n+(j-1)(alpha+1)+i;
		AppendTo[ rules, U[[j]][t+i-1] -> Z[[l]][t] ];
		AppendTo[ extEq, Z[[l]][t+1] == 
			If[ i==alpha+1, V[[j]][t], Z[[l+1]][t] ]
		],
	{j,m}, {i,alpha+1} ];
        {extEq /. rules /. rules1, rules }
] (* End Module *)

GeneralizedStateToExtendedState[ ___ ] := err /; Message[ GeneralizedStateToExtendedState::usage ]

Options[ GeneralizedStateToExtendedState ] = { NewVariables -> Automatic }



Lower[ stEq_, {X_, u_/;Head[u]=!=List}, t_, opts___] := 
	Lower[ stEq, {X, {u}}, t, opts ]

Lower[ stEq:{(_[t_+1]==_)..}, {X:{__}, U:{__}}, t_, opts___Rule] := 
	Module[ { lo, alpha, beta, met, n, m, Z, V, extEq, Hk, rules, transf, 
	transfinv, result, unabletointegrate },
	
	lo = LoweringOrder /. {opts} /. Options[Lower];
	met = Method /. {opts} /. Options[Lower];
	n = Length@X;
	alpha = Max[ Cases[ stEq, #[t+i_.]->i, -1]& /@ U];
	m = Length@U;
	beta = If[lo === Automatic, 0, alpha-lo ];

Which[ 
	Head[met] === List && Length[met] == n,
		transf = met,
	
	True,
		{extEq, rules} = GeneralizedStateToExtendedState[ stEq, {X, U}, t];
		Z = Array[ ToExpression["z" <> ToString[#] ]&, n + m(alpha+1) ];
		V = Array[ ToExpression["v" <> ToString[#] ]&, m ];
		Hk = SequenceH[ extEq, {Z, V}, t, alpha-beta+2 ];
    	transf = IntegratePfaff[
    		If[ lo === Automatic,
  				First@Select[ Reverse[Hk], Integrability, 1 ],
				Last[Hk]
    		]
    	];
		If[ FreeQ[ transf, De ], 
			transf = transf /. ((#[0] -> #[t]& /@ Z) /. Reverse /@ rules)
		]
]; (* End Which *)

	(* *************************** Main program ****************************** *)
result = Which[
	transf === {},
		{},
	!FreeQ[ transf, De],
		unabletointegrate = 1,
	True,
		transf = Complement[ transf, Flatten[ Cases[transf, #[_] ]&/@U ] ];
		transfinv = First@Solve[ OverTilde[#][t]& /@ X == transf,
			#[t]& /@ X ] /.
			c_Real?(Round[#] == #&) -> Round[c];
			(* Avaldan nihutamata v\[OTilde]rranditest x1[t],... *)
			Chop[
		Simplify[ First@ 
			Solve[ stEq /. (transfinv /. t -> t+1), 
				OverTilde[#][t+1]& /@ X 
			] (* End Solve *) /. transfinv 
		] (*End Simplify *), 
		epsilon ] (* End Chop*) /.
		c_Real?(Round[#] == #&) -> Round[c] /. 
		Rule -> Equal /. OverTilde[x_] -> x
];  (* End Which *)
result /; Not[ ValueQ[ unabletointegrate ] ]
] (* End Module*)

Lower[ ___ ] := err /; Message[ Lower::usage ]

Options[ Lower ] = { Method -> Automatic, LoweringOrder -> Automatic };



ClassicTransformability[ stEq_, {X_, u_/;Head[u]=!=List}, t_, opts___] := 
	ClassicTransformability[ stEq, {X, {u}}, t, opts ]
	
ClassicTransformability[ stEq:{(_[t_+1]==_)..}, {X_List, U_List}, t_, opts___] :=
Module[{n, alpha, m, Z, V, extEq },
	n = Length@X;
	alpha = Max[ Cases[ stEq, #[t+i_.]->i, -1]& /@ U];
	m = Length@U;
	met = Method /. {opts} /. Options[ ClassicTransformability ];
	Z = Array[ ToExpression["z" <> ToString[#] ]&, n + m(alpha+1) ];
    V = Array[ ToExpression["v" <> ToString[#] ]&, m ];
	extEq = First@GeneralizedStateToExtendedState[ stEq, {X, U}, t];
	Integrability[ Last[
		#[ extEq, {Z, V}, t, alpha+2 ]
	]] & [ If[ met === SubSpaceI, SequenceI, SequenceH ] ]
] (* End Module *)

ClassicTransformability[___] := err /; Message[ ClassicTransformability::usage ]

Options[ ClassicTransformability ] := { Method -> SubSpaceH }



ClassicTransform[ stEq_, {X_, u_/;Head[u]=!=List}, t_, opts___] :=
	ClassicTransform[ stEq, {X, {u}}, t, opts]

ClassicTransform[ stEq:{__}, {X:{__}, U:{__}}, t_, opts___Rule] := 
With[ { alpha = Max[ Cases[ stEq, #[t+i_.]->i, -1]& /@ U] },
	Lower[ stEq, {X, U}, t, opts, LoweringOrder -> alpha ]
]

ClassicTransform[ ___ ] := err /; Message[ ClassicTransform::usage ]

Options[ ClassicTransform ] = {Method -> Automatic}


On[General::spell1]

End[]

Protect[{
BackwardShift, BackwardShiftOperator, Basis, ClassicTransform,
ClassicTransformability, Coordinates, De, DifferentialForm,
Equivalence, ForwardShift, GeneralizedStateToExtendedState,
Integrability, IntegratePfaff, IntersectionSpace, IOToExtendedState,
IOToGeneralizedState, Irreducibility, Lower, LoweringOrder,
NewVariables, Realizability, Realization, ReducedDifferentialForm,
Reduction, SequenceH, SequenceI, ShowInfo, SimplifyBasis, Span,
SubmersionQ, SubSpaceH, SubSpaceI, ToDifferential, VerifySequence, 
Wedge, \[EmptySet]}
]

EndPackage[]