(* :Title: Utility functions for the package ExactModeling *) (* :Context: EMUtils` *) (* :Author: Martin Rickli Automatic Control Lab ETH Zurich, Physikstr. 3 CH-8092 Zurich Switzerland e-mail Rickli@aut.ee.ethz.ch *) (* :Summary: *) (* :Package Version: 1.0 *) (* :Mathematica Version: 2.2 *) (* :Copyright: Copyright 1993-96, Automatic Control Laboratory, ETH Zurich, Switzerland *) (* :History: 1.0 May 27 1996 (ri) final release *) (* :Keywords: Exact Modeling, Linear Systems, Interpolation, Identification, Polynomial Matrix *) (* :Sources: A.C. Antoulas & J.C. Willems: A behavioral approach to linear exact modeling, Transactions on Automatic Control, 38/12 p 1776-1802 (Dez. 1993) *) (* :Discussion: *) BeginPackage["EMUtils`"]; Off[General::spell] Off[General::spell1] (* --------------- >>> Usage Messages <<< --------------- *) (* usage messages for the exported functions and the context itself sorted alphabetically; Format: after '\n' 'Return' + 2 Spaces, no Space before '\n' *) EMUtils::usage = "EMUtils implements utility functions for use with package ExactModeling."; Options[EMUtils] = Sort[{ PrintLevel -> 3, PrintFlag -> False, Sorting -> MaxDifference, VectorNorm -> False, WorkingPrecision -> $MachinePrecision, }]; Assign::usage = "Assign[A, alist, row, col] overwrites elements in A at positions {row,col} with alist.\n Attributes[Assign] = {HoldFirst}.\n example:\n A = IdentityMatrix[8];\n alist = {{a, b, c}, {e, f, g}, {x, y, z}, {p, q, s}};\n row = {2, 4, 5, 7};\n col = {3, 6, 8};\n Assign[A, alist, row, col];"; CalculateCase::usage = "CalculateCase[{r,k},rowdeg] returns 'case'-number.\n case 0 : r = k = 0, \tzero vector\n case 1 : k = 0, \treal vector\n case 3 : rowdeg[[k]] - rowdeg[[r]] > 1\n case 5 : rowdeg[[r]] == rowdeg[[k]]\n case 7 : rowdeg[[k]] == rowdeg[[r]] + 1"; EMPrint::usage = "EMPrint[s] prints only on user request (PrintFlag->True)"; Freq2Data::usage = "Freq2Data[tf,lam] returns list of lists of data vectors w used in ProcessNewMeasurement.\n for each value in lam a matrix with as many rows as coloumns of tf is calculated.\n tf must be a pure function."; Householder::usage = "Householder[y] returns {h,a}, where h = I - 2(v.v')/(v'.v), householder-matrix a = Sqrt[y'.y], 2-norm of y such that h.y = [a 0 ... 0] Vector y must be real!"; Index::usage = "Index[expr] returns positions {r,k} of first Non-Zero element in expr. r is related to real part and k to imaginary. Vector expr is scaled first, so r > k."; LeadingCoefficientMatrix::usage = "LeadingCoefficientMatrix[m, var] returns the highest-row-degree coefficient matrix of polynomial matrix m.\n If this matrix is non-singular the matrix m is row reduced (see Kailath p 384)."; Linear::usage = "Sorting -> Linear sorts data in ascending order (first real second imaginary part)."; MaxDifference::usage = "Sorting -> MaxDifference sorts data such that difference of consecutive frequency values is as big as possible."; PrintInfo::usage = "PrintInfo[s] prints information if PrintLevel > 3"; RelatedTimeSeries::usage = "RelatedTimeSeries[p,lam,time] returns a related time series (eq. 4.1 in Ant_Paper)\n w_n,1 = (d/dtime - lam) ( p(time)Exp(lam time) )."; RowDeg::usage = "RowDeg[expr,var] returns the row degrees of the polynomial Matrix expr with indeterminate var."; RowReducedQ::usage = "RowReducedQ[m, var] returns True if polynomial matrix m(var) is row reduced, i.e. its leading coefficient matrix has full row rank."; ShowVariables::usage = "ShowVariables[context_String] prints all symbols defined in 'context' and their values."; SortFrequencies::usage = "SortFrequencies[vec, opts] Sorts vec such that adjacent elements differ 'a lot'. This may be used to determine sequence of calculations. The more consecutive frequencies differ the better numerical results can be expected."; SortIndex::usage = "SortIndex[vec] returns list of index needed to sort \"vec\" applying standard sorting procedure.\n sortedvec = vec[[indexlist]]"; SortPM::usage = "SortPM[th, var] returns {thnew, sortind} which reordered th in ascending row degrees."; SquareMatrixQ::usage = "SquareMatrixQ[m] is true if m is a square matrix."; ToNumber::usage = "ToNumber[num_, (prec), (opts)] sets small numbers to zero."; Assign::usage = "Assign[A, alist, row, col]"; MySingularValues::usage = "MySingularValues[m]"; PrintBugInfo::usage = "PrintBugInfo[strings, versionstring]"; PreSortIndex::usage = "PreSortIndex[err,rowdeg, cset]"; ToMatrixForm::usage = "ToMatrixForm[tf]"; TfData::usage = "TfData[fun, lam]"; Mean::usage = "Mean[list] returns the mean value of the elements in list."; NZeroVectorQ::usage = "NZeroVectorQ[vec, prec] returns True if elements of vec are smaller than 10^-prec.\n Default setting for prec is WorkingPrecision-WPOffset. See Options[EMUtils]"; NZeroQ::usage = "NZeroQ[n, prec]"; ZeroVectorQ::usage = "ZeroVectorQ[vec] returns True if vec === 0*vec."; Begin["`Private`"]; (* --------------- >>> Private Context <<< --------------- *) WPOffset = 5; (* Reduce WorkingPrecision when testing for zero *) EMUtils::nyimpl = "Not yet implemented case for `1`.\n\t`2`"; Attributes[ZeroQ] = {Listable}; ZeroQ[n_] := (n == 0); ZeroVectorQ[vec_] := ZeroQ[vec] /. List->And; (* Test if vec is identical to zero vector. Equal[] handles also the case of real numbers of the form 0. or 0.00 -> Zero *) Attributes[NZeroQ] = {Listable}; NZeroQ[n_, Infinity] := Together[n] === 0; NZeroQ[n_. sym_Symbol, prec_] := NZeroQ[n,prec]; NZeroQ[n_, prec_] := (Abs[N[n,prec] ] <= N[10^-(prec-5 (* WPOffset *) ),prec] ); NZeroQ[n_] := NZeroQ[n, (WorkingPrecision /. Options[EMUtils]) ]; NZeroVectorQ[vec_, prec_] := NZeroQ[vec,prec] /. List->And; NZeroVectorQ[vec_ ] := NZeroVectorQ[vec, (WorkingPrecision /. Options[EMUtils])]; (* handles also WorkingPrecision -> Infinity (<=> ZeroQ) /. List->And handles "all" structures => absolute test of each element. *) AdjustPrecision[num_, prec_] := N[Rationalize[num,0], prec]; AdjustPrecision::usage = "adjusts precision of 'num' to prec >= $MachinePrecision; prec = Infinity allowed"; ToNumber[num_, Infinity, opts___Rule] := Rationalize[Together @ ExpandAll[num], 0] (* works for numbers and vectors! *) ToNumber[num_, opts___Rule] := ToNumber[num, (WorkingPrecision /. {opts} /. Options[EMUtils]), opts] ToNumber[num_, precision_, opts___Rule] := (* converts num to an approximate number representation. precision = Infinity allowed. num may be a vector/matrix. values < 10^-(precision - WPOffset) set to zero (absolute test). If num contains symbolic expressions -> test every element. For vectors: VectorNorm->True uses norm of num (relative test). *) Module[{expr, prec, norm, dx, usenorm}, prec = Round[precision]; (* an Integer for use in Chop, N *) If[MemberQ[num, _Symbol, Infinity], usenorm = False, If[NZeroVectorQ[num,prec], Return[0 num] ]; usenorm = VectorNorm /. {opts} /. Options[EMUtils] ]; expr = N[ Together @ ExpandAll @ num, prec]; dx = 10^-(prec - WPOffset); If[usenorm && VectorQ[expr], norm = Sqrt @ Re[expr . Conjugate[expr] ]; dx = dx norm; (* -> rel. error bound *) ]; If[ Positive[dx], (* migth be zero if expr is zero or prec Infinity *) (* Chop small numbers; treats Re and Im seperately! *) expr = Chop[expr, dx] ]; expr ]; SquareMatrixQ[m_?MatrixQ] := SameQ @@ Dimensions[m]; (* copied from standard Package: LinearAlgebra`MatrixManipulation *) SquareMatrixQ[m_List] := False; Mean[list_List] := (Plus @@ list)/Length[list]; StandardDeviation[list_?VectorQ] := Sqrt[ Apply[Plus, (list - Mean[list])^2] / Length[list] ]; RowReducedQ[m_, var_Symbol] := !NZeroQ @ Det[ LeadingCoefficientMatrix[m,var] ]; Attributes[Assign] = { HoldFirst }; Assign[A_, alist_, row_, col_] := (* code from : Levent, lk3a@kelvin.seas.virginia.edu (via Mathgroup mailing list) Assign[A, alist, row, col] overwrites elements in A at positions {row, col} with alist. example: A = IdentityMatrix[8]; alist = {{a, b, c}, {e, f, g}, {x, y, z}, {p, q, s}}; row = {2, 4, 5, 7}; col = {3, 6, 8}; Assign[A, alist, row, col]; MatrixForm[A] *) Module[{blist, i, j, index}, index = Flatten[ Outer[ List, row, col ], 1 ]; blist = Flatten[alist]; Scan[ ({i, j} = index[[#]]; A[[i, j]] = blist[[#]])&, Range[Length[blist]] ] ]; (* Due to a bug in Mac-Mma 2.2 we have to redefine SingularValues. Suggestion from Dave Withoff, Wolfram Research (personal e-mail) *) MySingularValues[m_] := With[{result = SingularValues[m]}, If[Head[result] === List, result, N[SingularValues[SetPrecision[m, Precision[m] + 1]]]]] Householder[y_?VectorQ] /; !MemberQ[Head /@ y, Complex] := (* returns {h,a} h = I - 2(v.v')/(v'.v), householder-matrix a = Sqrt[y'.y], 2-norm of y Vector y must be real! such that h.y = [a 0 ... 0] h = h', h.h' = h'.h = I, h.h = I *) Module[{v=y,y1=First[y],a,h}, a = Sqrt[ Plus @@ (y y) ]; (* 2-norm of y *) a = - Sign[y1] a; v[[1]] = y1 - a; (* h = IdentityMatrix[Length[y] ] - 2 (Transpose[{v}].{v})/(Plus @@ (v v)); bug (?) in Dot: Precision lost if an element is zero! *) h = IdentityMatrix[Length[y] ] - 2 (Outer[Times,v,v]) / (Plus @@ (v v)); EMPrint["Householder {a,v} = ",{a,v}]; h = Together @ ExpandAll @ h; EMPrint[" h = ",h]; {h,a} ]; Householder[y_] := (Message[EMUtils::nyimpl,"Householder", "Argument must be a real valued vector."]; {0 y,0} (* Return (fictive) values {h,a} *) ); ShowVariables[context_String] := Module[{}, Print["Variables in context: "<>context]; Print["---------------------"]; Print[#,"\t",ToExpression[#] ]& /@ Names[context<>"*"]; ]; EMPrint[s___] := Print[s] /; (PrintFlag /. Options[EMUtils]); (* Printing should only occur on user request (PrintFlag->True) *) PrintInfo := If[(PrintLevel /. Options[EMUtils]) >= 3, Print ]; PrintBugInfo[st__String, ver_String] := Module[{}, Print[""]; Print["+--> Problem in ",First[{st}]," detected <--"]; Print["|"]; Print["| ",#]& /@ { "To help correct this problem we ask you to mail all", "calculations leading to this message to the author.", " E-mail: Rickli@aut.ee.ethz.ch", "Current version number: "<>ver, "System Information : "<>$System, "Mathematica Info : "<>$Version, "Thank you.",""}; If[Length[{st}] > 1, Print["+--> Additional information:\n|"]; Print["| ",#]& /@ Rest[{st}]; Print["|"]; ]; Print["+--<\n"]; ]; SortIndex[vec_] := (* 10. Feb. 1994, Ri Returns list of index needed to sort "vec" applying standard sorting procedure (ascending order!) sortedvec = vec[[indvec]] *) Module[{indvec}, indvec=Sort[ MapIndexed[List,vec] ]; Flatten @ Take[Transpose[indvec],{2}] ]; RowDeg[expr_,var_Symbol] := Module[{}, Flatten[ Map[ Max, Exponent[expr,var] ] ] ]; LeadingCoefficientMatrix[m_, var_Symbol] := Module[{rd}, rd = RowDeg[m,var]; MapIndexed[Coefficient[ ExpandAll[#1], var, rd[[First@#2]] ]&, m] ]; SortPM[th_, var_Symbol] := (* Returns {thnew, sortind} which reordered polynomial matrix th in ascending row degrees. *) Module[{thnew=th,rowdeg,sortind}, rowdeg = RowDeg[th,var]; sortind = SortIndex[rowdeg]; If[OrderedQ[rowdeg], Return[{th,sortind}] ]; thnew = th[[sortind]]; {thnew,sortind} ]; SortFrequencies[vec_List, opts___Rule] := (* Sort vec such that adjacent elements differ 'a lot'. This should be used to specify sequence of calculations. The more consecutive frequencies differ the better numerical results can be expected. Only useful in TF2Theta! *) Module[{len,v,sort}, len = Length[vec]; If[len < 2, Return[vec] ]; sort = Sorting /. {opts} /. Options[ExactModeling]; v = Sort[vec]; Switch[sort, None, v = vec, Linear, (* already sorted *), MaxDifference, If[EvenQ[len], v = Flatten @ Transpose[Partition[v,len/2] ], v = Flatten[ {Last[v], Transpose[Partition[v,(len-1)/2] ] }] ], _, v = vec ]; v ]; RelatedTimeSeries[p_,lam_, time_Symbol] := (* Returns vector p_n,1 of a related time series (eq. 4.1 in Ant_Paper) w_n,1 = (d/dt - lam) ( p(t)Exp(lam t) ) = dp/dt Exp(lam t). *) D[p,time]; Index[expr_?ZeroVectorQ] := {0, 0}; Index[expr_] := (* Returns index r and k: r refers to first nonzero element, k refers to first complex element. expr must have been scaled before (ScaleVector). Due to 'normalisation' index k is always greater than r! *) Module[{vec=expr,r=k=0}, r = Flatten @ Position[Re[vec], (p_/; p != 0)]; r = First[r]; (* r = 0 handled by testing ZeroVectorQ[vec] *) k = Flatten @ Position[Im[vec], (p_/; p != 0)]; k = Select[k,# > r&]; k = If[Length[k] == 0, 0, First[k] ]; {r,k} ]; CalculateCase::fatal = "Row degrees not ordered or invalid Indices r >= k > 0.\n rowdeg = `1`, {r,k} = `2`"; CalculateCase[{r_,k_},rowdeg_] /; (!OrderedQ[rowdeg] || (r >= k > 0) ):= (Message[CalculateCase::fatal,rowdeg,{r,k}]; Return[$Failed] ); CalculateCase[{0,0},rowdeg_] := (EMPrint["CalculateCase {r,k,case}: ",{0,0,0}];0); CalculateCase[{r_,0},rowdeg_] := (EMPrint["CalculateCase {r,k,case}: ",{r,0,1}];1); CalculateCase[{r_,k_},rowdeg_] := (* Assume th has ordered rowdegrees same case numbers as in Matlab files *) Module[{case}, Which[ (* r == 0 && k == 0,case = 0, k == 0, case = 1, *) rowdeg[[k]] - rowdeg[[r]] > 1, case = 3, rowdeg[[r]] == rowdeg[[k]], case = 5, rowdeg[[k]] == rowdeg[[r]] + 1, case = 7 ]; EMPrint["CalculateCase {r,k,case}: ",{r,k,case}]; case ]; PreSortIndex[err_?ZeroVectorQ, rowdeg_, cset_] := Range[Length[err] ]; PreSortIndex[err_,rowdeg_, cset_] := (* return index vector to reorder err and th, such that row r is in cset. affected rows have same row degree! err normalised and chopped (e.g. ToNumber). *) Module[{samedeg,r,ind=Range[Length[err]]}, (* try to reorder rows such that row r is in ControlSet *) r = Flatten @ Position[MapIndexed[#1=!=0&,err],True]; If[MemberQ[ControlSet,First[r]], Return[ind]]; samedeg = Select[r,rowdeg[[#]] == rowdeg[[First[r] ]]& ]; samedeg = Intersection[samedeg, cset]; r = First[r]; (* samedeg=list of row indices, index in cset error =!= 0 same row degree as row r r =!= samedeg (r NOT in cset) *) If[Length[samedeg] > 0, samedeg = First[samedeg]; ind = ind /. {r -> samedeg, samedeg -> r}; MyPrint["Exchange Rows for Preprocess: ",r,"<->",samedeg]; ]; ind ]; ToMatrixForm[tf_List] := (* return tf as a matrix *) Module[{}, If[!MatrixQ[tf], If[VectorQ[tf], (* assume single input *) Partition[tf,{1}], Partition[Flatten[{{tf}}], {1}] (* convert to single input! *) ], tf ] ]; ToMatrixForm[tf_] := ToMatrixForm[{tf}]; TfData[fun_Function, lam_List] := (* compute values of tf at lam. *) Module[{tf, data}, (* returns data as a List of Matrices *) tf = Function @@ { ToMatrixForm[Sequence @@ fun] }; data = (tf /@ lam) ]; Freq2Data[tf_Function, lam_List] := (* compute values of tf at lam. returns data as list of matrices (for each lam); nr of rows = nr of input in tf suitable for use in ProcessNewMeasurement. *) Module[{data,w,nin}, data = TfData[tf,lam]; nin = Dimensions[First @ data]; w = Transpose[Join[IdentityMatrix[nin[[2]] ],#] ]& /@ data ]; End[]; (* end the private context *) EndPackage[]; (* end the package context *) On[General::spell]; On[General::spell1];