wneuper/isa: doc-src/isac/jrocnik/Inverse_Z_Transform/Inverse_Z

     1 (* Title:  Test_Z_Transform

     2    Author: Jan Rocnik

     3    (c) copyright due to lincense terms.

     4 12345678901234567890123456789012345678901234567890123456789012345678901234567890

     5         10        20        30        40        50        60        70        80

     6 *)

     8 theory Inverse_Z_Transform imports Isac begin

    10 section {*trials towards Z transform *}

    11 text{*===============================*}

    12 subsection {*terms*}

    13 ML {*

    14 @{term "1 < || z ||"};

    15 @{term "z / (z - 1)"};

    16 @{term "-u -n - 1"};

    17 @{term "-u [-n - 1]"}; (*[ ] denotes lists !!!*)

    18 @{term "z /(z - 1) = -u [-n - 1]"};Isac

    19 @{term "1 < || z || ==> z / (z - 1) = -u [-n - 1]"};

    20 term2str @{term "1 < || z || ==> z / (z - 1) = -u [-n - 1]"};

    21 *}

    22 ML {*

    23 (*alpha -->  "</alpha>" *)

    25 @{term "\<alpha> "};

    26 @{term "\<delta> "};

    27 @{term "\<phi> "};

    28 @{term "\<rho> "};

    29 term2str @{term "\<rho> "};

    30 *}

    32 subsection {*rules*}

    33 (*axiomatization "z / (z - 1) = -u [-n - 1]" Illegal variable name: "z / (z - 1) = -u [-n - 1]" *)

    34 (*definition     "z / (z - 1) = -u [-n - 1]" Bad head of lhs: existing constant "op /"*)

    35 axiomatization where

    36   rule1: "1 = \<delta>[n]" and

    37   rule2: "|| z || > 1 ==> z / (z - 1) = u [n]" and

    38   rule3: "|| z || < 1 ==> z / (z - 1) = -u [-n - 1]" and

    39   rule4: "|| z || > || \<alpha> || ==> z / (z - \<alpha>) = \<alpha>^n * u [n]" and

    40   rule5: "|| z || < || \<alpha> || ==> z / (z - \<alpha>) = -(\<alpha>^n) * u [-n - 1]" and

    41   rule6: "|| z || > 1 ==> z/(z - 1)^2 = n * u [n]"

    42 ML {*

    43 @{thm rule1};

    44 @{thm rule2};

    45 @{thm rule3};

    46 @{thm rule4};

    47 *}

    49 subsection {*apply rules*}

    50 ML {*

    51 val inverse_Z = append_rls "inverse_Z" e_rls

    52   [ Thm  ("rule3",num_str @{thm rule3}),

    53     Thm  ("rule4",num_str @{thm rule4}),

    54     Thm  ("rule1",num_str @{thm rule1})

    55   ];

    57 val t = str2term "z / (z - 1) + z / (z - \<alpha>) + 1";

    58 val SOME (t', asm) = rewrite_set_ thy true inverse_Z t;

    59 term2str t' = "z / (z - ?\<delta> [?n]) + z / (z - \<alpha>) + ?\<delta> [?n]"; (*attention rule1 !!!*)

    60 *}

    61 ML {*

    62 val (thy, ro, er) = (@{theory}, tless_true, eval_rls);

    63 *}

    64 ML {*

    65 val SOME (t, asm1) = rewrite_ thy ro er true (num_str @{thm rule3}) t;

    66 term2str t = "- ?u [- ?n - 1] + z / (z - \<alpha>) + 1"; (*- real *)

    67 term2str t;

    68 *}

    69 ML {*

    70 val SOME (t, asm2) = rewrite_ thy ro er true (num_str @{thm rule4}) t;

    71 term2str t = "- ?u [- ?n - 1] + \<alpha> ^ ?n * ?u [?n] + 1"; (*- real *)

    72 term2str t;

    73 *}

    74 ML {*

    75 val SOME (t, asm3) = rewrite_ thy ro er true (num_str @{thm rule1}) t;

    76 term2str t = "- ?u [- ?n - 1] + \<alpha> ^ ?n * ?u [?n] + ?\<delta> [?n]"; (*- real *)

    77 term2str t;

    78 *}

    79 ML {*

    80 terms2str (asm1 @ asm2 @ asm3);

    81 *}

    83 section {*Prepare steps in CTP-based programming language*}

    84 text{*===================================================*}

    85 subsection {*prepare expression*}

    86 ML {*

    87 val ctxt = ProofContext.init_global @{theory};

    88 val ctxt = declare_constraints' [@{term "z::real"}] ctxt;

    90 val SOME fun1 = parseNEW ctxt "X z = 3 / (z - 1/4 + -1/8 * z ^ -1)"; term2str fun1;

    91 val SOME fun1' = parseNEW ctxt "X z = 3 / (z - 1/4 + -1/8 * (1/z))"; term2str fun1';

    92 *}

    94 axiomatization where

    95   ruleZY: "(X z = a / b) = (X' z = a / (z * b))"

    97 ML {*

    98 val (thy, ro, er) = (@{theory}, tless_true, eval_rls);

    99 val SOME (fun2, asm1) = rewrite_ thy ro er true  @{thm ruleZY} fun1; term2str fun2;

   100 val SOME (fun2', asm1) = rewrite_ thy ro er true  @{thm ruleZY} fun1'; term2str fun2';

   102 val SOME (fun3,_) = rewrite_set_ @{theory Isac} false norm_Rational fun2;

   103 term2str fun3; (*fails on x^(-1) TODO*)

   104 val SOME (fun3',_) = rewrite_set_ @{theory Isac} false norm_Rational fun2';

   105 term2str fun3'; (*OK*)

   107 val (_, expr) = HOLogic.dest_eq fun3'; term2str expr;

   108 *}

   110 subsection {*solve equation*}

   111 ML {*(*from test/Tools/isac/Minisubpbl/100-init-rootpbl.sml*)

   112 "----------- Minisubplb/100-init-rootp (*OK*)bl.sml ---------------------";

   113 val denominator = parseNEW ctxt "z^2 - 1/4*z - 1/8 = 0";

   114 val fmz = ["equality (z^2 - 1/4*z - 1/8 = (0::real))", "solveFor z","solutions L"];

   115 val (dI',pI',mI') =("Isac", ["univariate","equation"], ["no_met"]);

   116 (*                           ^^^^^^^^^^^^^^^^^^^^^^ TODO: ISAC determines type of eq*)

   117 *}

   118 ML {*

   119 val (p,_,f,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))];

   120 val (p,_,f,nxt,_,pt) = me nxt p [] pt;

   121 (*[

   122 (([], Frm), solve (z ^ 2 - 1 / 4 * z - 1 / 8 = 0, z)),

   123 (([1], Frm), z ^ 2 - 1 / 4 * z - 1 / 8 = 0),              bad rewrite order

   124 (([1], Res), -1 / 8 + z ^ 2 + -1 / 4 * z = 0),            bad pattern

   125 (([2], Pbl), solve (-1 / 8 + z ^ 2 + -1 / 4 * z = 0, z)),

   126 (([2,1], Pbl), solve (-1 / 8 + z ^ 2 + -1 / 4 * z = 0, z)),

   127 (([2,1,1], Pbl), solve (-1 / 8 + z ^ 2 + -1 / 4 * z = 0, z)),

   128 (([2,1,1,1], Frm), -1 / 8 + z ^ 2 + -1 / 4 * z = 0)]

   129 *)

   130 *}

   131 ML {*

   132 val denominator = parseNEW ctxt "-1/8 + -1/4*z + z^2 = 0";

   133 (*ergebnis: [gleichung, was tun?, lösung]*)

   134 val fmz = ["equality (-1/8 + -1/4*z + z^2 = (0::real))", "solveFor z","solutions L"];

   135 (*liste der theoreme die zum lösen benötigt werden, aus isac, keine spezielle methode (no met)*)

   136 val (dI',pI',mI') =

   137   ("Isac", ["pqFormula","degree_2","polynomial","univariate","equation"], ["no_met"]);

   138 (*schritte abarbeiten*)

   139 val (p,_,f,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))];

   140 val (p,_,f,nxt,_,pt) = me nxt p [] pt;

   141 val (p,_,f,nxt,_,pt) = me nxt p [] pt;

   142 val (p,_,f,nxt,_,pt) = me nxt p [] pt;

   143 val (p,_,f,nxt,_,pt) = me nxt p [] pt; (*val nxt = ("Empty_Tac", ...): tac'_*)

   144 show_pt pt;

   145 *}

   147 subsection {*partial fraction decomposition*}

   148 subsubsection {*solution of the equation*}

   149 ML {*

   150 val SOME solutions = parseNEW ctxt "[z=1/2, z=-1/4]";

   151 term2str solutions;

   152 atomty solutions;

   153 *}

   155 subsubsection {*get solutions out of list*}

   156 text {*in isac's CTP-based programming language: $let s_1 = NTH 1 solutions; s_2 = NTH 2...$*}

   157 ML {*

   158 val Const ("List.list.Cons", _) $ s_1 $ (Const ("List.list.Cons", _) $

   159       s_2 $ Const ("List.list.Nil", _)) = solutions;

   160 term2str s_1;

   161 term2str s_2;

   162 *}

   164 ML {* (*Solutions as Denominator --> Denominator1 = z - Zeropoint1, Denominator2 = z-Zeropoint2,...*)

   165 val xx = HOLogic.dest_eq s_1;

   166 val s_1' = HOLogic.mk_binop "Groups.minus_class.minus" xx;

   167 val xx = HOLogic.dest_eq s_2;

   168 val s_2' = HOLogic.mk_binop "Groups.minus_class.minus" xx;

   169 term2str s_1';

   170 term2str s_2';

   171 *}

   173 subsubsection {*build expression*}

   174 text {*in isac's CTP-based programming language: $let s_1 = Take numerator / (s_1 * s_2)$*}

   175 ML {*

   176 (*The Main Denominator is the multiplikation of the partial fraction denominators*)

   177 val denominator' = HOLogic.mk_binop "Groups.times_class.times" (s_1', s_2') ;

   178 val SOME numerator = parseNEW ctxt "3::real";

   180 val expr' = HOLogic.mk_binop "Rings.inverse_class.divide" (numerator, denominator');

   181 term2str expr';

   182 *}

   184 subsubsection {*Ansatz - create partial fractions out of our expression*}

   186 axiomatization where

   187   ansatz2: "n / (a*b) = A/a + B/(b::real)" and

   188   multiply_eq2: "(n / (a*b) = A/a + B/b) = (a*b*(n  / (a*b)) = a*b*(A/a + B/b))"

   190 ML {*

   191 (*we use our ansatz2 to rewrite our expression and get an equilation with our expression on the left and the partial fractions of it on the right side*)

   192 val SOME (t1,_) = rewrite_ @{theory Isac} e_rew_ord e_rls false @{thm ansatz2} expr';

   193 term2str t1;

   194 atomty t1;

   195 val eq1 = HOLogic.mk_eq (expr', t1);

   196 term2str eq1;

   197 *}

   198 ML {*

   199 (*eliminate the demoninators by multiplying the left and the right side with the main denominator*)

   200 val SOME (eq2,_) = rewrite_ @{theory Isac} e_rew_ord e_rls false @{thm multiply_eq2} eq1;

   201 term2str eq2;

   202 *}

   203 ML {*

   204 (*simplificatoin*)

   205 val SOME (eq3,_) = rewrite_set_ @{theory Isac} false norm_Rational eq2;

   206 term2str eq3; (*?A ?B not simplified*)

   207 *}

   208 ML {*

   209 val SOME fract1 =

   210   parseNEW ctxt "(z - 1 / 2) * (z - -1 / 4) * (A / (z - 1 / 2) + B / (z - -1 / 4))"; (*A B !*)

   211 val SOME (fract2,_) = rewrite_set_ @{theory Isac} false norm_Rational fract1;

   212 term2str fract2 = "(A + -2 * B + 4 * A * z + 4 * B * z) / 4";

   213 (*term2str fract2 = "A * (1 / 4 + z) + B * (-1 / 2 + z)" would be more traditional*)

   214 *}

   215 ML {*

   216 val (numerator, denominator) = HOLogic.dest_eq eq3;

   217 val eq3' = HOLogic.mk_eq (numerator, fract1); (*A B !*)

   218 term2str eq3';

   219 *}

   220 ML {* (*MANDATORY: otherwise 3 = 0*)

   221 val SOME (eq3'' ,_) = rewrite_set_ @{theory Isac} false norm_Rational eq3';

   222 term2str eq3'';

   223 *}

   225 subsubsection {*get first koeffizient*}

   227 ML {*

   228 (*substitude z with the first zeropoint to get A*)

   229 val SOME (eq4_1,_) = rewrite_terms_ @{theory Isac} e_rew_ord e_rls [s_1] eq3'';

   230 term2str eq4_1;

   231 *}

   232 ML {*

   233 val SOME (eq4_2,_) = rewrite_set_ @{theory Isac} false norm_Rational eq4_1;

   234 term2str eq4_2;

   235 *}

   236 ML {*

   237 val fmz = ["equality (3 = 3 * A / (4::real))", "solveFor A","solutions L"];

   238 val (dI',pI',mI') =("Isac", ["univariate","equation"], ["no_met"]);

   240 *}

   241 ML {*

   242 (*solve the simple linear equilation for A TODO: return eq, not list of eq*)

   243 val (p,_,fa,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))];

   244 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   245 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   246 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   247 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   248 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   249 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   250 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   251 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   252 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   253 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   254 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   255 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   256 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   257 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   258 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   259 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   260 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   261 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   262 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   263 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   264 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   265 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   266 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   267 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   268 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   269 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   270 *}

   271 ML {*

   272 val (p,_,fa,nxt,_,pt) = me nxt p [] pt;

   273 f2str fa;

   274 *}

   276 subsubsection {*get second koeffizient*}

   278 ML {*

   279 (*substitude z with the second zeropoint to get B*)

   280 val SOME (eq4b_1,_) = rewrite_terms_ @{theory Isac} e_rew_ord e_rls [s_2] eq3'';

   281 term2str eq4b_1;

   282 *}

   284 ML {*

   285 val SOME (eq4b_2,_) = rewrite_set_ @{theory Isac} false norm_Rational eq4b_1;

   286 term2str eq4b_2;

   287 *}

   289 ML {*

   290 (*solve the simple linear equilation for B TODO: return eq, not list of eq*)

   291 val fmz = ["equality (3 = -3 * B / (4::real))", "solveFor B","solutions L"];

   292 val (dI',pI',mI') =("Isac", ["univariate","equation"], ["no_met"]);

   293 val (p,_,fb,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))];

   294 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   295 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   296 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   297 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   298 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   299 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   300 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   301 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   302 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   303 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   304 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   305 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   306 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   307 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   308 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   309 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   310 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   311 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   312 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   313 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   314 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   315 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   316 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   317 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   318 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   319 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   320 val (p,_,fb,nxt,_,pt) = me nxt p [] pt;

   321 f2str fb;

   322 *}

   324 ML {* (*check koeffizients*)

   325 if f2str fa = "[A = 4]" then () else error "part.fract. eq4_1";

   326 if f2str fb = "[B = -4]" then () else error "part.fract. eq4_1";

   327 *}

   329 subsubsection {*substitute expression with solutions*}

   330 ML {*

   331 *}

   333 section {*Implement the Specification and the Method*}

   334 text{*==============================================*}

   335 subsection{*Define the Specification*}

   336 ML {*

   337 val thy = @{theory};

   338 *}

   339 ML {*

   340 store_pbt

   341  (prep_pbt thy "pbl_SP" [] e_pblID

   342  (["SignalProcessing"], [], e_rls, NONE, []));

   343 store_pbt

   344  (prep_pbt thy "pbl_SP_Ztrans" [] e_pblID

   345  (["Z_Transform","SignalProcessing"], [], e_rls, NONE, []));

   346 store_pbt

   347  (prep_pbt thy "pbl_SP_Ztrans_inv" [] e_pblID

   348  (["inverse", "Z_Transform", "SignalProcessing"],

   349   [("#Given" ,["equality X_eq"]),

   350    ("#Find"  ,["equality n_eq"])

   351   ],

   352   append_rls "e_rls" e_rls [(*for preds in where_*)], NONE,

   353   [["TODO: path to method"]]));

   355 show_ptyps();

   356 get_pbt ["inverse","Z_Transform","SignalProcessing"];

   357 *}

   359 subsection{*Define the (Dummy-)Method*}

   360 subsection {*Define Name and Signature for the Method*}

   361 consts

   362   InverseZTransform :: "[bool, bool] => bool"

   363     ("((Script InverseZTransform (_ =))// (_))" 9)

   365 ML {*

   366 store_met

   367  (prep_met thy "met_SP" [] e_metID

   368  (["SignalProcessing"], [],

   369    {rew_ord'="tless_true", rls'= e_rls, calc = [], srls = e_rls, prls = e_rls,

   370     crls = e_rls, nrls = e_rls}, "empty_script"));

   371 store_met

   372  (prep_met thy "met_SP_Ztrans" [] e_metID

   373  (["SignalProcessing", "Z_Transform"], [],

   374    {rew_ord'="tless_true", rls'= e_rls, calc = [], srls = e_rls, prls = e_rls,

   375     crls = e_rls, nrls = e_rls}, "empty_script"));

   376 *}

   377 ML {*

   378 store_met

   379  (prep_met thy "met_SP_Ztrans_inv" [] e_metID

   380  (["SignalProcessing", "Z_Transform", "inverse"], [],

   381    {rew_ord'="tless_true", rls'= e_rls, calc = [], srls = e_rls, prls = e_rls,

   382     crls = e_rls, nrls = e_rls},

   383   "empty_script"

   384  ));

   385 *}

   386 ML {*(*

   387 store_met

   388  (prep_met thy "met_SP_Ztrans_inv" [] e_metID

   389  (["SignalProcessing", "Z_Transform", "inverse"], [],

   390    {rew_ord'="tless_true", rls'= e_rls, calc = [], srls = e_rls, prls = e_rls,

   391     crls = e_rls, nrls = e_rls},

   392   "Script InverseZTransform (Xeq::bool) =" ^

   393   " (let X = Take Xeq;" ^

   394   "      X = Rewrite ruleZY False X" ^

   395   "  in X)"

   396  ));

   397 *)*}

   398 ML {*

   399 show_mets();

   400 get_met ["SignalProcessing","Z_Transform","inverse"];

   401 *}

   404 section {*Program in CTP-based language*}

   405 text{*=================================*}

   406 subsection {*Stepwise extend Program*}

   407 ML {*

   408 val str =

   409 "Script InverseZTransform (Xeq::bool) =" ^

   410 " Xeq";

   411 *}

   412 ML {*

   413 val str =

   414 "Script InverseZTransform (Xeq::bool) =" ^

   415 " (let X = Take Xeq;" ^

   416 "      X = Rewrite ruleZY False X" ^

   417 "  in X)";

   418 *}

   419 ML {*

   420 val thy = @{theory};

   421 val sc = ((inst_abs thy) o term_of o the o (parse thy)) str;

   422 *}

   423 ML {*

   424 term2str sc;

   425 atomty sc

   426 *}

   428 subsection {*Stepwise Execute the Program*}

   437 section {*Write Tests for Crucial Details*}

   438 text{*===================================*}

   439 ML {*

   441 *}

   443 section {*Integrate Program into Knowledge*}

   444 ML {*

   446 *}

   448 end

author	Jan Rocnik <jan.rocnik@student.tugraz.at>
	Thu, 08 Sep 2011 23:17:35 +0200
branch	decompose-isar
changeset 42252	e633bb41ea42
permissions	-rwxr-xr-x