source: git/Singular/LIB/control.lib @ 0266ac

version="$Id: control.lib,v 1.33 2005-05-23 15:47:35 levandov Exp $";
category="System and Control Theory";
info="
LIBRARY:  control.lib     Algebraic analysis tools for System and Control Theory

AUTHORS:  Oleksandr Iena       yena@mathematik.uni-kl.de
@*        Markus Becker        mbecker@mathematik.uni-kl.de
@*        Viktor Levandovskyy  levandov@mathematik.uni-kl.de

SUPPORT: Forschungsschwerpunkt 'Mathematik und Praxis' (Project of Dr. E. Zerz
and V. Levandovskyy), Uni Kaiserslautern

MAIN PROCEDURES:
  control(R);            analysis of controllability-related properties of R (using Ext modules)
  controlDim(R);           analysis of controllability-related properties of R (using dimension)
  autonom(R);            analysis of autonomy-related properties of R (using Ext modules)
  autonomDim(R);           analysis of autonomy-related properties of R (using dimension)

COMPONENT PROCEDURES:
  leftKernel(R);         a left kernel of R
  rightKernel(R);        a right kernel of R
  leftInverse(R);        a left inverse of R
  rightInverse(R);       a right inverse of R
  smith(M);              a Smith form of a module M
  colrank(M);            a column rank of M as of matrix
  genericity(M);         analysis of the genericity of parameters
  canonize(L);           Groebnerification for modules in the output of control or autonomy procs
  iostruct(R);           computes an IO-structure of behavior given by a module R
  findTorsion(R, I);     generators of the submodule of a module R, annihilated by the ideal I

AUXILIARY PROCEDURES:
  controlExample(s);     set up an example from the mini database inside of the library
  view();                well-formatted output of lists, modules and matrices
";

LIB "homolog.lib";
LIB "poly.lib";
LIB "primdec.lib";
LIB "matrix.lib";

//---------------------------------------------------------------
static proc Opt_Our()
"USAGE: Opt_Our();
RETURN: intvec, where previous options are stored
PURPOSE: save previous options and set customized options
"
{
  intvec v;
  v=option(get);
  option(redSB);
  option(redTail);
  return (v);
}

//-------------------------------------------------------------------------

static proc space(int n)
"USAGE:spase(n); n is an integer (number of needed spaces)
RETURN:  string consisting of n spaces
NOTE:  the procedure is used in the procedure 'view' to have a better formatted output
"{
  int i;
  string s="";
  for(i=1;i<=n;i++)
  {
    s=s+" ";
  };
  return(s);
};
//-----------------------------------------------------------------------------
proc view(M)
"USAGE:  view(M);   M is of any type
RETURN:  no return value
PURPOSE:  procedure for (well-) formatted output of modules, matrices, lists of modules, matrices; shows everything even if entries are long
NOTE:  in case of other types( not 'module', 'matrix', 'list') works just as standard 'print' procedure
EXAMPLE:  example view; shows an example
"{
  // to be replaced with something more feasible
  if ( (typeof(M)=="module")||(typeof(M)=="matrix") )
  {
  int @R=nrows(M);
  int @C=ncols(M);
  int i;
  int j;
  list MaxLength=list();
  int Size=0;
  int max;
  string s;

  for(i=1;i<=@C;i++)
  {
    max=0;

    for(j=1;j<=@R;j++)
    {
      Size=size( string( M[j,i] ) );
      if( Size>max )
      {
        max=Size;
      };
    };
    MaxLength[i] = max;
  };

  for(i=1;i<=@R;i++)
  {
    s="";
    for(j=1;j<@C;j++)
    {
      s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) ) +",";
    };

    s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) );

    if (i!=@R)
    {
      s=s+",";
    };
    print(s);
  };

  return();
  };

  if(typeof(M)=="list")
  {
    int sz=size(M);
    int i;
    for(i=1;i<=sz;i++)
    {
      view(M[i]);
      print("");
    };

    return();
  };
  print(M);
  return();
}
example
{"EXAMPLE:";echo = 2;
  ring r;
  list L;
  matrix M[1][3] = x2+x,y3-y,z5-4z+7;
  L[1] = "a matrix:"; 
  L[2] = M;
  L[3] = "an ideal:";
  L[4] = ideal(M);
  view(L);
};
//--------------------------------------------------------------------------
proc rightKernel(matrix M)
"USAGE:  rightKernel(M);   M a matrix
RETURN:  module
PURPOSE: computes the right kernel of matrix M (a module of all elements v such that Mv=0)
EXAMPLE: example rightKernel; shows an example
"{
  return(modulo(M,std(0)));
}
example
{"EXAMPLE:";echo = 2;
  ring r = 0,(x,y,z),dp;
  matrix M[1][3] = x,y,z;
  print(M);
  matrix R = rightKernel(M);
  print(R);
  // check:
  print(M*R);
};
//-------------------------------------------------------------------------
proc leftKernel(matrix M)
"USAGE:  leftKernel(M);   M a matrix
RETURN:  module
PURPOSE: computes left kernel of matrix M (a module of all elements v such that vM=0)
EXAMPLE:  example leftKernel; shows an example
"
{
  return( transpose( modulo( transpose(M),std(0) ) ) );
}
example
{"EXAMPLE:";echo = 2;
  ring r= 0,(x,y,z),dp;
  matrix M[3][1] = x,y,z;
  print(M);
  matrix L = leftKernel(M);
  print(L);
  // check:
  print(L*M);
};
//------------------------------------------------------------------------
proc leftInverse(module M)
"USAGE:  leftInverse(M);  M a module
RETURN:  module
PURPOSE: computes such a matrix L, that LM = Id;
EXAMPLE: example leftInverse; shows an example
NOTE: exists only in the case when M is free submodule
"
{
  // it works also for the NC case;
  int NCols = ncols(M);
  module Id = freemodule(NCols);
  module N  = transpose(M);
  intvec old_opt=Opt_Our();
  Id = std(Id);
  matrix T;
  // check the correctness (Id \subseteq M)
  // via dimension: dim (M) = -1!
  int d = dim_Our(N);
  if (d != -1)
  {
    // No left inverse exists
    return(matrix(0));
  }
  matrix T2 = lift(N, Id);
  T2 =  transpose(T2);
  option(set,old_opt);  // set the options back
  return(T2);
}
example
{
  "EXAMPLE:";echo =2;
  // a trivial example:
  ring r = 0,(x,z),dp;
  matrix M[2][1] = 1,x2z;
  print(M);
  print( leftInverse(M) );
  kill r;
  // derived from the example TwoPendula:
  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  matrix U[3][1];
  U[1,1]=(-L2)*Dt^4+(g)*Dt^2;
  U[2,1]=(-L1)*Dt^4+(g)*Dt^2;
  U[3,1]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2);
  module M = module(U);
  module L = leftInverse(M);
  print(L);
  // check
  print(L*M);
};
//-----------------------------------------------------------------------
proc rightInverse(module R)
"USAGE:  rightInverse(M);  M a module
RETURN:  module
PURPOSE: computes such a matrix L, that ML = Id
EXAMPLE: example rightInverse; shows an example
NOTE: exists only in the case when M is free submodule
"
{
  return(transpose(leftInverse(transpose(R))));
}
example
{ "EXAMPLE:";echo =2;
  // a trivial example:
  ring r = 0,(x,z),dp;
  matrix M[1][2] = 1,x2+z;
  print(M);
  print( rightInverse(M) );
  kill r;
  // derived from the TwoPendula example:
  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  matrix U[1][3];
  U[1,1]=(-L2)*Dt^4+(g)*Dt^2;
  U[1,2]=(-L1)*Dt^4+(g)*Dt^2;
  U[1,3]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2);
  module M = module(U);
  module L = rightInverse(M);
  print(L);
  // check
  print(M*L);
};
//-----------------------------------------------------------------------
static proc dim_Our(module R)
{
  int d;
  if (attrib(R,"isSB")<>1)
  {
    R = std(R);
  }
  d = dim(R);
  return(d);
}
//-----------------------------------------------------------------------
static proc Ann_Our(module R)
{
  return(Ann(R));
}
//-----------------------------------------------------------------------
static proc Ext_Our(int i, module R, list #)
{
  // mimicking 'Ext_R' from homolog.lib
  int ExtraArg = ( size(#)>0 );
  if (ExtraArg)
  {
    return( Ext_R(i,R,#[1]) );
  }
  else
  {
    return( Ext_R(i,R) );
  }
}
//------------------------------------------------------------------------
static proc is_zero_Our
{
  //just a copy of 'is_zero' from "poly.lib" patched with GKdim
  int d = dim_Our(std(#[1]));
  int a = ( d==-1 );
  if( size(#) >1 ) { list L=a,d; return(L); }
  return(a);
  //  return( is_zero(R) ) ;
};
//------------------------------------------------------------------------
static proc control_output(int i, int NVars, module R, module Ext_1, list Gen)
"USAGE:  control_output(i, NVars, R, Ext_1),
PURPOSE: where
@*         i is integer (number of first nonzero Ext or a number of variables in a basering + 1 in case that all the Exts are zero),
@*           NVars:  integer, number of variables in a base ring,
@*           R:  module R (cokernel representation),
@*           Ext_1:  module, the first Ext(its cokernel representation)
RETURN:  list with all the contollability properties of the system which is to be returned in 'control' procedure
NOTE:  this procedure is used in 'control' procedure
"{
  // TODO: NVars to be replaced with the global hom. dimension of basering!!!
  // Is not clear what to do with gl.dim of qrings
  string DofS = "dimension of the system:";
  string Fn   = "number of first nonzero Ext:";
  string Gen_mes = "Parameter constellations which might lead to a non-controllable system:";

  module RK = rightKernel(R);
  int d=dim_Our(std(transpose(R)));

  if (i==1)
  {
    return(
            list ( Fn,
                   i,
                  "not controllable , image representation for controllable part:",
                    RK,
                  "kernel representation for controllable part:",
                   leftKernel( RK ),
                  "obstruction to controllability",
                   Ext_1,
                  "annihilator of torsion module (of obstruction to controllability)",
                   Ann_Our(Ext_1),
                   DofS,
                   d
                 )
          );
  };

  if(i>NVars)
  {
    return( list(  Fn,
                   -1,
                  "strongly controllable(flat), image representation:",
                   RK,
                  "left inverse to image representation:",
                   leftInverse(RK),
                   DofS,
                   d,
                   Gen_mes,
                   Gen)
          );
  };

  //
  //now i<=NVars
  //

  if( (i==2) )
  {
    return( list( Fn,
                  i,
                 "controllable, not reflexive, image representation:",
                  RK,
                  DofS,
                  d,
                  Gen_mes,
                  Gen)
          );
  };

  if( (i>=3) )
  {
    return( list ( Fn,
                   i,
                  "reflexive, not strongly controllable, image representation:",
                   RK,
                   DofS,
                   d,
                      Gen_mes,
                   Gen)
          );
  };
};
//-------------------------------------------------------------------------

proc control(module R)
"USAGE:  control(R);  R a module (R is the matrix of the system of equations to be investigated)
RETURN:  list
PURPOSE: compute the list of all the properties concerning controllability of the system (behavior), represented by the matrix R
EXAMPLE:  example control; shows an example
"
{
  int i;
  int NVars=nvars(basering);
  // TODO: NVars to be replaced with the global hom. dimension of basering!!!
  int ExtIsZero;
  intvec v=Opt_Our();
  module R_std=std(R);
  module Ext_1 = std(Ext_Our(1,R_std));

  ExtIsZero=is_zero_Our(Ext_1);
  i=1;
  while( (ExtIsZero) && (i<=NVars) )
  {
    i++;
    ExtIsZero = is_zero_Our( Ext_Our(i,R_std) );
  };
  matrix T=lift(R,R_std);
  list l=genericity(T);
  option(set,v);

  return( control_output( i, NVars, R, Ext_1, l ) );
}
example
{"EXAMPLE:";echo = 2;
  // a WindTunnel example
  ring A = (0,a, omega, zeta, k),(D1, delta),dp;
  module R;
  R = [D1+a, -k*a*delta, 0, 0],
      [0, D1, -1, 0],
      [0, omega^2, D1+2*zeta*omega, -omega^2];
  R=transpose(R);
  view(R);
  view(control(R));
};
//--------------------------------------------------------------------------
proc controlDim(module R)
"USAGE:  controlDim(R); R a module (R is the matrix of the system of equations to be investigated)
RETURN: list
PURPOSE: computes list of all the properties concerning controllability of the system (behavior), represented by the  matrix R
EXAMPLE:  example controlDim; shows an example
NOTE: this procedure is analogous to 'control' but uses dimension calculations.This approach works for full row rank matrices only.
"
{
  if( nrows(R) != colrank(transpose(R)) )
  {
    return ("controlDim cannot be applied, since R does not have full row rank");
  }
  intvec v=Opt_Our();
  module R_std=std(R);
  int d=dim_Our(R_std);
  int NVars=nvars(basering);
  int i=NVars-d;
  module Ext_1=std(Ext_Our(1,R_std));
  matrix T=lift(R,R_std);
  list l=genericity(T);
  option(set, v);
  return( control_output( i, NVars, R, Ext_1, l));
}
example
{"EXAMPLE:";echo = 2;
  //a WindTunnel example
  ring A = (0,a, omega, zeta, k),(D1, delta),dp;
  module R;
  R = [D1+a, -k*a*delta, 0, 0],
      [0, D1, -1, 0],
      [0, omega^2, D1+2*zeta*omega, -omega^2];
  R=transpose(R);
  view(R);
  view(controlDim(R));
};
//------------------------------------------------------------------------
proc colrank(module M)
"USAGE: colrank(M); M a matrix/module
RETURN: int
PURPOSE: compute the column rank of M as of matrix
NOTE: this procedure uses Bareiss algorithm 
"{
  // NOte continued:
  // which might not terminate in some cases
  module M_red  = bareiss(M)[1];
  int NCols_red = ncols(M_red);
  return (NCols_red);
}
example
{"EXAMPLE: ";echo = 2;
  // de Rham complex
  ring r=0,(D(1..3)),dp;
  module R;
  R=[0,-D(3),D(2)],
    [D(3),0,-D(1)],
    [-D(2),D(1),0];
  R=transpose(R);
  colrank(R);
};

//------------------------------------------------------------------------
static proc autonom_output( int i, int NVars, module RC, int R_rank )
"USAGE:  proc autonom_output(i, NVars, RC, R_rank)
           i:  integer, number of first nonzero Ext or
               just number of variables in a base ring + 1 in case that all the Exts are zero
           NVars:  integer, number of variables in a base ring
           RC: module, kernel-representation of controllable part of the system
           R_rank: integer, column rank of the representation matrix
PURPOSE: compute all the autonomy properties of the system which is to be returned in 'autonom' procedure
RETURN:  list
NOTE:  this procedure is used in 'autonom' procedure
"
{
  int d=NVars-i;//that is the dimension of the system
  string DofS="dimension of the system:";
  string Fn = "number of first nonzero Ext:";
  if(i==0)
  {
    return( list(  Fn,
                   i,
                  "not autonomous",
                   "kernel representation for controllable part",
                   RC,
                   "column rank of the matrix",
                   R_rank,
                   DofS,
                   d )
          );
  };

  if( i>NVars )
  {
    return( list( Fn,
                  -1,
                  "trivial",
                  DofS,
                  d )
          );
  };

  //
  //now i<=NVars
  //


  if( i==1 )
  // in case that NVars==1 there is no sense to consider the notion
  // of strongly autonomous behavior, because it does not imply
  // that system is overdetermined in this case
  {
    return( list ( Fn,
                   i,
                  "autonomous, not overdetermined",
                   DofS,
                   d )
          );
  };

  if( i==NVars )
  {
    return( list(  Fn,
                   i,
                  "strongly autonomous(fin. dimensional),in particular overdetermined",
                   DofS,
                   d)
          );
  };

  if( i<NVars )
  {
    return( list ( Fn,
                   i,
                  "overdetermined, not strongly autonomous",
                   DofS,
                   d)
          );
  };
};
//--------------------------------------------------------------------------
proc autonomDim(module R)
"USAGE:  autonomDim(R);   R a module (R is a matrix of the system of equations which is to be investigated)
RETURN: list
PURPOSE: computes the list of all the properties concerning autonomy of the system (behavior), represented by the matrix R
NOTE:  this procedure is analogous to 'autonom' but uses dimension calculations
EXAMPLE:  example autonomDim; shows an example
"
{
  int d;
  int NVars = nvars(basering);
  module RT = transpose(R);
  module RC;  //for computation of controllable part if if exists
  int R_rank = ncols(R);
  d     = dim_Our( std(RT) );  //this is the dimension of the system
  int i = NVars-d;  //First non-zero Ext
  if( d==0 )
    {
      RC=leftKernel(rightKernel(R));
      R_rank=colrank(R);
    }
  return( autonom_output(i,NVars,RC,R_rank) );
}
example
{"EXAMPLE:"; echo = 2;
  // Cauchy1 example
  ring r=0,(s1,s2,s3,s4),dp;
  module R= [s1,-s2],
            [s2, s1],
            [s3,-s4],
            [s4, s3];
  R=transpose(R);
  view( R );
  view( autonomDim(R) );
};
//----------------------------------------------------------
proc autonom(module R)
"USAGE:  autonom(R);   R a module (R is a matrix of the system of equations which is to be investigated)
RETURN:  list
PURPOSE: find all the properties concerning autonomy of the system (behavior) represented by the  matrix R
EXAMPLE: example autonom; shows an example
"
{
  int NVars=nvars(basering);
  int ExtIsZero;
  module RT=transpose(R);
  module RC;
  int R_rank=ncols(R);
  ExtIsZero=is_zero_Our(Ext_Our(0,RT));
  int i=0;
  while( (ExtIsZero)&&(i<=NVars) )
  {
    i++;
    ExtIsZero = is_zero_Our(Ext_Our(i,RT));
  };
  if (i==0)
  {
    RC=leftKernel(rightKernel(R));
    R_rank=colrank(R);
  }
  return(autonom_output(i,NVars,RC,R_rank));
}
example
{"EXAMPLE:"; echo = 2;
  // Cauchy
  ring r=0,(s1,s2,s3,s4),dp;
  module R= [s1,-s2],
            [s2, s1],
            [s3,-s4],
            [s4, s3];
  R=transpose(R);
  view( R );
  view( autonom(R) );
};


//----------------------------------------------------------
proc genericity(matrix M)
"USAGE:  genericity(M); M is a matrix/module
RETURN:  list (of strings)
PURPOSE: determine parametric expressions which have been assumed to be non-zero in the process of computing the Groebner basis
NOTE: we strongly recommend to switch on the redSB and redTail options;
@*    the procedure is effective with the lift procedure for modules with parameters
EXAMPLE:  example genericity; shows an example
"
{
  // returns "-", if there are no parameters!
  if (npars(basering)==0)
  {
    return("-");
  }
  list RT = evas_genericity(M); // list of strings
  if ((size(RT)==1) && (RT[1] == ""))
  {
    return("-");
  }
  return(RT);
}
example
{  // TwoPendula
  "EXAMPLE:"; echo = 2;
  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  module RR =
    [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2],
    [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
    [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];
  module R = transpose(RR);
  module SR = std(R);
  matrix T = lift(R,SR);
  genericity(T);
  //-- The result might be different when computing reduced bases:
  matrix T2;
  option(redSB);
  option(redTail);
  module SR2 = std(R);
  T2 =  lift(R,SR2);
  genericity(T2);
}
//---------------------------------------------------------------
static proc victors_genericity(matrix M)
{
  // returns "-", if there are no parameters!
  if (npars(basering)==0)
  {
    return("-");
  }
  int plevel = printlevel-voice+2;
  // M is a matrix over a ring with params and vars;
  ideal I = ideal(M); // a list of entries
  I = simplify(I,2); // delete 0's
  // decompose every coeff in every poly
  int i;
  int s = size(I);
  ideal NM;
  poly p;
  number num;
  int  cl=1;
  intvec ZeroVec; ZeroVec[nvars(basering)] = 0;
  intvec W;
  ideal Numero, Denomiro;
  int cNu=0; int cDe=0;
  for (i=1; i<=s; i++)
  {
    // remove contents and add them as polys
    p   = I[i];
    W   = leadexp(p);
    if (W == ZeroVec) // i.e. just a coef
    {
      num    = denominator(leadcoef(p)); // from poly.lib
      NM[cl] = numerator(leadcoef(p));
      dbprint(p,"numerator:");
      dbprint(p, string(NM[cl]));
      cNu++; Numero[cNu]= NM[cl];
      cl++;
      NM[cl] = num; // denominator
      dbprint(p,"denominator:");
      dbprint(p, string(NM[cl]));
      cDe++; Denomiro[cDe]= NM[cl];
      cl++;
      p = p - lead(p); // for the next cycle
    }
    if ( p!= 0)
    {
      num = content(p);
      p   = p/num;
      NM[cl] = denominator(num);
      dbprint(p,"content denominator:");
      dbprint(p, string(NM[cl]));
      cNu++; Numero[cNu]= NM[cl];
      cl++;
      NM[cl] = numerator(num);
      dbprint(p,"content numerator:");
      dbprint(p, string(NM[cl]));
      cDe++; Denomiro[cDe]= NM[cl];
      cl++;
    }
    // it seems that the next elements will not have real influence
    while( p != 0)
    {
      NM[cl] = leadcoef(p); // should be all integer, i.e. non-rational
      dbprint(p,"coef:");
      dbprint(p, string(NM[cl]));
      cl++;
      p = p - lead(p);
    }
  }
  NM = simplify(NM,4); // delete identical
  string newvars = parstr(basering);
  def save = basering;
  string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
  execute(NewRing);
  // get params as variables
  // create a list of non-monomials
  ideal @L;
  ideal F;
  ideal NM = imap(save,NM);
  NM = simplify(NM,8); //delete multiples
  poly p,q;
  cl = 1;
  int j, cf;
  for(i=1; i<=size(NM);i++)
  {
    p = NM[i] - lead(NM[i]);
    if (p!=0)
    {
      //      L[cl] = p;
      F  = factorize(NM[i],1); //non-constant factors only
      cf = 1;
      // factorize every polynomial
      // throw away every monomial from factorization (also constants from above ring)
      for (j=1; j<=size(F);j++)
      {
        q = F[j]-lead(F[j]);
        if (q!=0)
        {
          @L[cl] = F[j];
          cl++;
        }
      }
    }
  }
  // return the result [in string-format]
  @L = simplify(@L,2+4+8); // skip zeroes, doubled and entries, diff. by a constant
  list SL;
  for (j=1; j<=size(@L);j++)
  {
    SL[j] = string(@L[j]);
  }
  setring save;
  return(SL);
}
//---------------------------------------------------------------
static proc evas_genericity(matrix M)
{
  // called from the main genericity proc
  ideal I = ideal(M);
  I = simplify(I,2+4);
  int s = size(I);
  ideal Den;
  poly p;
  int i;
  for (i=1; i<=s; i++)
  {
    p = I[i];
    while (p !=0)
    {
      Den = Den, denominator(leadcoef(p));
      p   = p-lead(p);
    }
  }
  Den = simplify(Den,2+4);
  string newvars = parstr(basering);
  def save = basering;
  string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
  execute(NewRing);
  ideal F;
  ideal Den = imap(save,Den);
  Den = simplify(Den,2);
  int s1 = size(Den);
  for (i=1; i<=s1; i++)
  {
    if (Den[i] !=1)
    {
      F= F, factorize(Den[i],1);
    }
  }
  F = simplify(F, 2+4+8);
  ideal @L = F;
  list SL;
  int c,j;
  string Mono;
  c = 1;
  for (j=1; j<=size(@L);j++)
  {
    if (leadcoef(@L[j]) <0)
    {
      @L[j] = -1*@L[j];
    }
    if ( (@L[j] - lead(@L[j]))==0 ) //@L[j] is a monomial
    {
      Mono = Mono + string(@L[j])+ ","; // concatenation
    }
    else
    {
      c++;
      SL[c] = string(@L[j]);
    }
  }
  if (Mono!="")
  {
    Mono = Mono[1..size(Mono)-1]; // delete the last semicolon
  }
  SL[1] = Mono;
  setring save;
  return(SL);
}

//---------------------------------------------------------------
proc canonize(list L)
"USAGE:  canonize(L); L a list
RETURN:  list
PURPOSE: modules in the list are canonized by computing their reduced minimal (= unique up to constant factor w.r.t. the given ordering) Groebner bases
ASSUME:  L is the output of control/autonomy procedures
EXAMPLE:  example canonize; shows an example
"
{
  list M = L;
  intvec v=Opt_Our();
  int s = size(L);
  int i;
  for (i=2; i<=s; i=i+2)
  {
    if (typeof(M[i])=="module")
    {
      M[i] = std(M[i]);
      //      M[i] = prune(M[i]); // mimimal embedding: no need yet
      //      M[i] = std(M[i]);
    }
  }
  option(set, v); //set old values back
  return(M);
}
example
{  // TwoPendula with L1=L2=L
  "EXAMPLE:"; echo = 2;
  ring r=(0,m1,m2,M,g,L),Dt,dp;
  module RR =
    [m1*L*Dt^2, m2*L*Dt^2, -1, (M+m1+m2)*Dt^2],
    [m1*L^2*Dt^2-m1*L*g, 0, 0, m1*L*Dt^2],
    [0, m2*L^2*Dt^2-m2*L*g, 0, m2*L*Dt^2];
  module R = transpose(RR);
  list C = control(R);
  list CC = canonize(C);
  view(CC);
}

//----------------------------------------------------------------

static proc elementof (int i, intvec v)
{
  int b=0;
  for(int j=1;j<=nrows(v);j++)
    {
      if(v[j]==i)
        {
          b=1;
          return (b);
        }
    }
  return (b);
}
//-----------------------------------------------------------------
proc iostruct(module R)
"USAGE: iostruct( R ); R a module
RETURN:  list L with entries: string s, intvec v, module P and module Q
PURPOSE:  if R is the kernel-representation-matrix of some system, then we output a input-ouput representation Py=Qu of the system, the components that have been chosen as outputs(intvec v) and a comment s
NOTE:  the procedure uses Bareiss algorithm 
EXAMPLE:  example iostruct; shows an example
"
{
  // NOTE cont'd
  //which might not terminate in some cases
  list L = bareiss(R);
  int R_rank = ncols(L[1]);
  int NCols=ncols(R);
  intvec v=L[2];
  int temp;
  int NRows=nrows(v);
  int i,j;
  int b=1;
  module P;
  module Q;
  int n=0;

  while(b==1)               //sort v through bubblesort
    {
      b=0;
      for(i=1;i<NRows;i++)
        {
          if(v[i]>v[i+1])
          {
            temp=v[i];
            v[i]=v[i+1];
            v[i+1]=temp;
            b=1;
          }
        }
    }
  P=R[v];                     //generate P
  for(i=1;i<=NCols;i++)       //generate Q
    {
      if(elementof(i,v)==1)
        {
          i++;
          continue;
        }
      Q=Q,R[i];
    }
  Q=simplify(Q,2);
  string s="The following components have been chosen as outputs: ";
  return (list(s,v,P,Q));
}
example
{"EXAMPLE:";echo = 2;
 //Example Antenna
 ring r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), (c,dp);
 module RR;
 RR =
   [Dt, -K1, 0, 0, 0, 0, 0, 0, 0],
   [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta],
   [0, 0, Dt, -K1, 0, 0, 0, 0, 0],
   [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta],
   [0, 0, 0, 0, Dt, -K1, 0, 0, 0],
   [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta];
 module R = transpose(RR);
 view(iostruct(R));
};

//---------------------------------------------------------------
static proc smdeg(matrix N)
"USAGE: smdeg( N ); N a matrix
RETURN:  intvec
PURPOSE: returns an intvec of length 2 with the index of an element of N with smallest degree
"
{
  int n = nrows(N);
  int m = ncols(N);
  int d,d_temp;
  intvec v;
  int i,j;            // counter

  if (N==0)
  {
    v = 1,1;
    return(v);
  }

  for (i=1; i<=n; i++)
// hier wird ein Element ausgewaehlt(!=0) und mit dessen Grad gestartet
  {
    for (j=1; j<=m; j++)
    {
      if( deg(N[i,j])!=-1 )
      {
        d=deg(N[i,j]);
        break;
      }
    }
    if (d != -1)
    {
      break;
    }
  }
  for(i=1; i<=n; i++)
  {
    for(j=1; j<=m; j++)
    {
      d_temp = deg(N[i,j]);
      if ( (d_temp < d) && (N[i,j]!=0) )
      {
        d=d_temp;
      }
    }
  }
  for (i=1; i<=n; i++)
  {
    for (j=1; j<=m;j++)
    {
      if ( (deg(N[i,j]) == d) && (N[i,j]!=0) )
      {
        v = i,j;
        return(v);
      }
    }
  }
}
//---------------------------------------------------------------
static proc NoNon0Pol(vector v)
"USAGE: NoNon0Pol(v), v a vector
RETURN:  int
PURPOSE: returns 1, if there is only one non-zero element in v and 0 else
"{
  int i,j;
  int n = nrows(v);
  for( j=1; j<=n; j++)
  {
    if (v[j] != 0)
    {
      i++;
    }
  }
  if ( i!=1 )
  {
    i=0;
  }
  return(i);
}
//---------------------------------------------------------------
static proc extgcd_Our(poly p, poly q)
{
  ideal J;   //for extgcd-computations
  matrix T; //----------"------------
  list L;
  // the extgcd-command has a bug in versions before 2-0-7
  if ( system("version")<=2006 )
  {
    J = p,q; // J = N[k-1,k-1],N[k,k]; //J is of type ideal
    L[1] = liftstd(J,T);  //T is of type matrix
    if(J[1]==p) //this is just for the case the SINGULAR swaps the
    //      two elements due to ordering
    {
      L[2] = T[1,1];
      L[3] = T[2,1];
    }
    else
    {
      L[2] = T[2,1];
      L[3] = T[1,1];
    }
  }
  else
  {
    L=extgcd(p,q);
    //    L=extgcd(N[k-1,k-1],N[k,k]);
    //one can use this line if extgcd-bug is fixed
  }
  return(L);
}
static proc normalize_Our(matrix N, matrix Q)
"USAGE: normalize_Our(N,Q), N, Q are two matrices
PURPOSE: normalizes N and divides the columns of Q through the leading coefficients of the columns of N
RETURN: normalized matrix N and altered Q(according to the scheme mentioned in purpose). If number of columns of N and Q do not coincide, N and Q are returned unchanged
NOTE: number of columns of N and Q must coincide.
"
{
  if(ncols(N) != ncols(Q))
    {
      return (N,Q);
    }
  module M = module(N);
  module S = module(Q);
  int NCols = ncols(N);
  number n;
  for(int i=1;i<=NCols;i++)
    {
      n = leadcoef(M[i]);
      if( n != 0 )
        {
          M[i]=M[i]/n;
          S[i]=S[i]/n;
        }
     }
   N = matrix(M);
   Q = matrix(S);
   return (N,Q);
}

//---------------------------------------------------------------
proc smith( module M )
"USAGE: smith(M); M a module/matrix
PURPOSE: computes the Smith normal form of a matrix
RETURN: a list of length 4 with the following entries:
@*      [1]: the Smith normal form S of M,
@*      [2]: the rank of M,
@*      [3]: a unimodular matrix U,
@*      [4]: a unimodular matrix V,
such that U*M*V=S. An warning is returned when no Smith form exists.
NOTE: The Smith form only exists over PIDs (principal ideal domains). Use global ordering for computations!
"
{
  if (nvars(basering)>1) //if more than one variable, return empty list
  {
    string s="The Smith-Form only exists for principal ideal domains";
    return (s);
  }
  matrix N = matrix(M);         //Typecasting
  int n = nrows(N);
  int m = ncols(N);
  matrix P = unitmat(n);       //left transformation matrix
  matrix Q = unitmat(m);       //right transformation matrix
  int k, i, j, deg_temp;
  poly tmp;
  vector v;
  list L;                      //for extgcd-computation
  intmat f[1][n];              //to save degrees
  matrix lambda[1][n];         //to save leadcoefficients
  intmat g[1][m];              //to save degrees
  matrix mu[1][m];             //to save leadcoefficients
  int ii;                       //counter

  while ((k!=n) && (k!=m) )
  {
    k++;
    while ((k<=n) && (k<=m))  //outer while-loop for column-operations
    {
      while(k<=m )        //inner while-loop for row-operations
      {
        if( (n>m) && (k < n) && (k<m))
        {
          if( simplify((ideal(submat(N,k+1..n,k+1..m))),2)== 0)
          {
            return(N,k-1,P,Q);
          }
        }
        i,j = smdeg(submat(N,k..n,k..m)); //choose smallest degree in the remaining submatrix
        i=i+(k-1);                        //indices adjusted to the whole matrix
        j=j+(k-1);
        if(i!=k)                    //take the element with smallest degree in the first position
        {
          N=permrow(N,i,k);
          P=permrow(P,i,k);
        }
        if(j!=k)
        {
          N=permcol(N,j,k);
          Q=permcol(Q,j,k);
        }
        if(NoNon0Pol(N[k])==1)
        {
          break;
        }
        tmp=leadcoef(N[k,k]);
        deg_temp=ord(N[k,k]);             //ord outputs the leading degree of N[k,k]
        for(ii=k+1;ii<=n;ii++)
        {
          lambda[1,ii]=leadcoef(N[ii,k])/tmp;
          f[1,ii]=deg(N[ii,k])-deg_temp;
        }
        for(ii=k+1;ii<=n;ii++)
        {
          N = addrow(N,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
          P = addrow(P,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
          N,Q=normalize_Our(N,Q);
        }
      }
      if (k>n)
      {
        break;
      }
      if(NoNon0Pol(transpose(N)[k])==1)
      {
        break;
      }
      tmp=leadcoef(N[k,k]);
      deg_temp=ord(N[k,k]); //ord outputs the leading degree of N[k][k]

      for(ii=k+1;ii<=m;ii++)
      {
        mu[1,ii]=leadcoef(N[k,ii])/tmp;
        g[1,ii]=deg(N[k,ii])-deg_temp;
      }
      for(ii=k+1;ii<=m;ii++)
      {
        N=addcol(N,k,-mu[1,ii]*var(1)^g[1,ii],ii);
        Q=addcol(Q,k,-mu[1,ii]*var(1)^g[1,ii],ii);
        N,Q=normalize_Our(N,Q);
      }
    }
    if( (k!=1) && (k<n) && (k<m) )
    {
      L = extgcd_Our(N[k-1,k-1],N[k,k]);
      if ( N[k-1,k-1]!=L[1] )  //means that N[k-1,k-1] is not a divisor of N[k,k]
      {
        N=addrow(N,k-1,L[2],k);
        P=addrow(P,k-1,L[2],k);
        N,Q=normalize_Our(N,Q);

        N=addcol(N,k,-L[3],k-1);
        Q=addcol(Q,k,-L[3],k-1);
        N,Q=normalize_Our(N,Q);
        k=k-2;
      }
    }
  }
  if( (k<=n) && (k<=m) )
  {
    if( N[k,k]==0)
    {
      return(N,k-1,P,Q);
    }
  }
  return(N,k,P,Q);
}
example
{
  "EXAMPLE:";echo = 2;
  option(redSB);
  option(redTail);
  ring r   = 0,x,dp;
  module M = [x2,x,3x3-4], [2x2-1,4x,5x2], [2x5,3x,4x];
  print(M);
  list P = smith(M);
  print(P[1]);
  matrix N = matrix(M);
  matrix B = P[3]*N*P[4];
  print(B);
}
// see what happens when the matrix is already in Smith-Form
//  module M = [x,0,0],[0,x2,0],[0,0,x3];
//  list L = smith(M);
// print(L[1]);
//matrix N=matrix(M);
//matrix B=L[3]*N*L[4];
//print(B);
//---------------------------------------------------------------
static proc list_tex(L, string name,link l,int nr_loop)
"USAGE: list_tex(L,name,l), where L is a list, name a string, l a link
         writes the content of list L in a tex-file 'name'
RETURN: nothing
"
{
  if(typeof(L)!="list")  //in case L is not a list
  {
    texobj(name,L);
  }
  if(size(L)==0)
  {
  }
  else
  {
    string t;
    for (int i=1;i<=size(L);i++)
    {
      while(1)
      {
        if(typeof(L[i])=="string")  //Fehler hier fuer normalen output->nur wenn string in liste dann verbatim
        {
          t=L[i];
          if(nr_loop==1)
          {
            write(l,"\\begin\{center\}");
            write(l,"\\begin\{verbatim\}");
          }
          write(l,t);
          if(nr_loop==0)
          {
            write(l,"\\par");
          }
          if(nr_loop==1)
          {
            write(l,"\\end\{verbatim\}");
            write(l,"\\end\{center\}");
          }
          break;
        }
        if(typeof(L[i])=="module")
        {
          texobj(name,matrix(L[i]));
          break;
        }
        if(typeof(L[i])=="list")
        {
          list_tex(L[i],name,l,1);
          break;
        }
        write(l,"\\begin\{center\}");
        texobj(name,L[i]);
        write(l,"\\end\{center\}");
        write(l,"\\par");
        break;
      }
    }
  }
}
example
{
  "EXAMPLE:";echo = 2;
}
//---------------------------------------------------------------
proc verbatim_tex(string s, link l)
"USAGE: verbatim_tex(s,l), where s is a string and l a link
PURPOSE: writes the content of s in verbatim-environment in the file
         specified by link
RETURN: nothing
"
{
  write(l,"\\begin{verbatim}");
  write(l,s);
  write(l,"\\end{verbatim}");
  write(l,"\\par");
}
example
{
  "EXAMPLE:";echo = 2;
}
//---------------------------------------------------------------
proc findTorsion(module R, ideal TAnn)
"USAGE:  findTorsion(R, I);   R an ideal/matrix/module, I an ideal
RETURN:  module
PURPOSE: computes the Groebner basis of the submodule of R, annihilated by I
NOTE: especially helpful, when I is the annihilator of the t(R) - the torsion submodule of R. In this case, the result is the explicit presentation of t(R) as
the submodule of R
EXAMPLE: example findTorsion; shows an example
"
{
  // motivation: let R be a module,
  // TAnn is the annihilator of t(R)\subset R
  // compute the generators of t(R) explicitly
  ideal AS = TAnn;
  module S = R;
  if (attrib(S,"isSB")<>1)
  {
    S = std(S);
  }
  if (attrib(AS,"isSB")<>1)
  {
    AS = std(AS);
  }
  int nc  = ncols(S);
  module To = quotient(S,AS);
  To = std(NF(To,S));
  return(To);
}
example
{
  "EXAMPLE:";echo = 2;
  // Flexible Rod
  ring A = 0,(D1, D2), (c,dp);
  module R= [D1, -D1*D2, -1], [2*D1*D2, -D1-D1*D2^2, 0];
  module RR = transpose(R);
  list L = control(RR);
  // here, we have the annihilator:
  ideal LAnn = D1; // = L[10]
  module Tr  = findTorsion(RR,LAnn);
  print(RR);  // the module itself
  print(Tr); // generators of the torsion submodule
}


proc controlExample(string s)
"USAGE:  controlExample(s);   s a string
RETURN:  ring
PURPOSE: set up an example from the mini database by initalizing a ring and a module in a ring
NOTE: in order to see the list of available examples, execute @code{controlExample(\"show\");}
@* To use ab example, one has to do the following. Suppose one calls the ring, where the example will be activated, A. Then, by executing
@*  @code{def A = controlExample(\"Antenna\");} and @code{setring A;},
@* A will become a basering from the example \"Antenna\" with
the predefined system module R (transposed).
After that one can just execute @code{control(R);} respectively
@code{autonom(R);} to perform the control resp. autonomy analysis of R.
EXAMPLE: example controlExample; shows an example
"{
  list E, S, D; // E=official name, S=synonym, D=description
  E[1] = "Cauchy1";  S[1] = "cauchy1";  D[1] = "1-dimensional Cauchy equation";
  E[2] = "Cauchy2";  S[2] = "cauchy2";  D[2] = "2-dimensional Cauchy equation";
  E[3] = "Control1"; S[3] = "control1"; D[3] = "example of a simple noncontrollable system";
  E[4] = "Control2"; S[4] = "control2"; D[4] = "example of a simple controllable system";
  E[5] = "Antenna";  S[5] = "antenna";  D[5] = "antenna";
  E[6] = "Einstein"; S[6] = "einstein"; D[6] = "Einstein equations in vacuum";
  E[7] = "FlexibleRod"; S[7] = "flexible rod"; D[7] = "flexible rod";
  E[8] = "TwoPendula";  S[8] = "two pendula";  D[8] = "two pendula mounted on a cart";
  E[9] = "WindTunnel";  S[9] = "wind tunnel";D[9] = "wind tunnel";
  E[10] = "Zerz1";      S[10] = "zerz1"; D[10] = "example from the lecture of Eva Zerz";
  // all the examples so far
  int i;
  if ( (s=="show") || (s=="Show") )
  {
    print("The list of examples:");
    for (i=1; i<=size(E); i++)
    {
      printf("name: %s,  desc: %s", E[i],D[i]);
    }
    return();
  }
  string t;
  for (i=1; i<=size(E); i++)
  {
    if ( (s==E[i]) || (s==S[i]) )
    {
      t = "def @A = ex"+E[i]+"();";
      execute(t);
      return(@A);
    }
  }
  "No example found";
  return();
}
example
{
  "EXAMPLE:";echo = 2;
  controlExample("show");   // let us see all available examples:
  def B = controlExample("TwoPendula"); // let us set up a particular example
  setring B;
  print(R);
}

//----------------------------------------------------------
//
//Some example rings with defined systems
//----------------------------------------------------------
//autonomy:
//
//----------------------------------------------------------
static proc exCauchy1()
{
  ring @r=0,(s1,s2),dp;
  module R= [s1,-s2],
            [s2, s1];
  R=transpose(R);
  export R;
  return(@r);
};
//----------------------------------------------------------
static proc exCauchy2()
{
  ring @r=0,(s1,s2,s3,s4),dp;
  module R= [s1,-s2],
            [s2, s1],
            [s3,-s4],
            [s4, s3];
  R=transpose(R);
  export R;
  return(@r);
};
//----------------------------------------------------------
static proc exZerz1()
{
  ring @r=0,(d1,d2),dp;
  module R=[d1^2-d2],
           [d2^2-1];
  R=transpose(R);
  export R;
  return(@r);
};
//----------------------------------------------------------
//control
//----------------------------------------------------------
static proc exControl1()
{
  ring @r=0,(s1,s2,s3),dp;
  module R=[0,-s3,s2],
           [s3,0,-s1];
  R=transpose(R);
  export R;
  return(@r);
};
//----------------------------------------------------------
static proc exControl2()
{
  ring @r=0,(s1,s2,s3),dp;
  module R=[0,-s3,s2],
           [s3,0,-s1],
           [-s2,s1,0];
  R=transpose(R);
  export R;
  return(@r);
};
//----------------------------------------------------------
static proc exAntenna()
{
  ring @r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), dp;
  module R;
  R = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0],
      [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta],
      [0, 0, Dt, -K1, 0, 0, 0, 0, 0],
      [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta],
      [0, 0, 0, 0, Dt, -K1, 0, 0, 0],
      [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta];

  R=transpose(R);
  export R;
  return(@r);
};

//----------------------------------------------------------

static proc exEinstein()
{
  ring @r = 0,(D(1..4)),dp;
  module R =
  [D(2)^2+D(3)^2-D(4)^2, D(1)^2, D(1)^2, -D(1)^2, -2*D(1)*D(2), 0, 0, -2*D(1)*D(3), 0, 2*D(1)*D(4)],
  [D(2)^2, D(1)^2+D(3)^2-D(4)^2, D(2)^2, -D(2)^2, -2*D(1)*D(2), -2*D(2)*D(3), 0, 0, 2*D(2)*D(4), 0],
  [D(3)^2, D(3)^2, D(1)^2+D(2)^2-D(4)^2, -D(3)^2, 0, -2*D(2)*D(3), 2*D(3)*D(4), -2*D(1)*D(3), 0, 0],
  [D(4)^2, D(4)^2, D(4)^2, D(1)^2+D(2)^2+D(3)^2, 0, 0, -2*D(3)*D(4), 0, -2*D(2)*D(4), -2*D(1)*D(4)],
  [0, 0, D(1)*D(2), -D(1)*D(2), D(3)^2-D(4)^2, -D(1)*D(3), 0, -D(2)*D(3), D(1)*D(4), D(2)*D(4)],
  [D(2)*D(3), 0, 0, -D(2)*D(3),-D(1)*D(3), D(1)^2-D(4)^2, D(2)*D(4), -D(1)*D(2), D(3)*D(4), 0],
  [D(3)*D(4), D(3)*D(4), 0, 0, 0, -D(2)*D(4), D(1)^2+D(2)^2, -D(1)*D(4), -D(2)*D(3), -D(1)*D(3)],
  [0, D(1)*D(3), 0, -D(1)*D(3), -D(2)*D(3), -D(1)*D(2), D(1)*D(4), D(2)^2-D(4)^2, 0, D(3)*D(4)],
  [D(2)*D(4), 0, D(2)*D(4), 0, -D(1)*D(4), -D(3)*D(4), -D(2)*D(3), 0, D(1)^2+D(3)^2, -D(1)*D(2)],
  [0, D(1)*D(4), D(1)*D(4), 0, -D(2)*D(4), 0, -D(1)*D(3), -D(3)*D(4), -D(1)*D(2), D(2)^2+D(3)^2];

  R=transpose(R);
  export R;
  return(@r);
};

//----------------------------------------------------------
static proc exFlexibleRod()
{
  ring @r = 0,(D1, delta), dp;
  module R;
  R = [D1, -D1*delta, -1], [2*D1*delta, -D1-D1*delta^2, 0];

  R=transpose(R);
  export R;
  return(@r);
};

//----------------------------------------------------------
static proc exTwoPendula()
{
  ring @r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  module R = [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2],
             [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
             [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];

  R=transpose(R);
  export R;
  return(@r);
};
//----------------------------------------------------------
static proc exWindTunnel()
{
  ring @r = (0,a, omega, zeta, k),(D1, delta),dp;
  module R = [D1+a, -k*a*delta, 0, 0],
             [0, D1, -1, 0],
             [0, omega^2, D1+2*zeta*omega, -omega^2];

  R=transpose(R);
  export R;
  return(@r);
};