test/Tools/isac/Knowledge/partial_fractions.sml
author Walther Neuper <neuper@ist.tugraz.at>
Thu, 05 Apr 2012 16:53:43 +0200
changeset 42405 f813ece49902
parent 42386 3aff35f94465
child 42413 a8303098408a
permissions -rwxr-xr-x
renamed ID "inverse" to "Inverse" (Isabelle2002 --> 2011)

reason: "inverse" is occupied by Const ("Rings.inverse_class.inverse", ...
     1 (* Title:  partial fraction decomposition
     2    Author: Jan Rocnik
     3    (c) due to copyright terms
     4 *)
     5 
     6 "--------------------------------------------------------";
     7 "table of contents --------------------------------------";
     8 "--------------------------------------------------------";
     9 "----------- why helpless here ? ------------------------";
    10 "----------- why not nxt = Model_Problem here ? ---------";
    11 "----------- fun factors_from_solution ------------------";
    12 "----------- Logic.unvarify_global ----------------------";
    13 "----------- eval_drop_questionmarks --------------------";
    14 "----------- = me for met_partial_fraction --------------";
    15 "----------- progr.vers.2: check erls for multiply_ansatz";
    16 "----------- progr.vers.2: improve program --------------";
    17 "--------------------------------------------------------";
    18 "--------------------------------------------------------";
    19 "--------------------------------------------------------";
    20 
    21 
    22 "----------- why helpless here ? ------------------------";
    23 "----------- why helpless here ? ------------------------";
    24 "----------- why helpless here ? ------------------------";
    25 val fmz = ["filterExpression (X z = 3 / (z - 1/4 + -1/8 * (1/(z::real))))", 
    26   "stepResponse (x[n::real]::bool)"];
    27 val (dI,pI,mI) = ("Isac", ["Inverse", "Z_Transform", "SignalProcessing"], 
    28   ["SignalProcessing","Z_Transform","Inverse"]);
    29 val (p,_,f,nxt,_,pt)  = CalcTreeTEST [(fmz, (dI,pI,mI))]; 
    30 val (p,_,f,nxt,_,pt) = me nxt p [] pt; "Add_Given";
    31 val (p,_,f,nxt,_,pt) = me nxt p [] pt; "Add_Find";
    32 val (p,_,f,nxt,_,pt) = me nxt p [] pt; "Specify_Theory";
    33 val (p,_,f,nxt,_,pt) = me nxt p [] pt; "Specify_Problem";
    34 val (p,_,f,nxt,_,pt) = me nxt p [] pt; "Specify_Method";
    35 val (p,_,f,nxt,_,pt) = me nxt p [] pt; "nxt = Apply_Method";
    36 val (p,_,f,nxt,_,pt) = me nxt p [] pt; "nxt = Rewrite (ruleZY, Inverse_Z_Transform.ruleZY) --> X z = 3 / (z - 1 / 4 + -1 / 8 * (1 / z))";
    37 val (p''',_,f,nxt''',_,pt''') = me nxt p [] pt; "nxt = Rewrite_Set norm_Rational --> X' z = 3 / (z * (z - 1 / 4 + -1 / 8 * (1 / z)))";
    38 "~~~~~ fun me, args:"; val ((_,tac), (p:pos'), _, (pt:ptree)) = (nxt, p, [], pt);
    39 val ("ok", (_, _, ptp)) = locatetac tac (pt,p)
    40 val (pt, p) = ptp;
    41 "~~~~~ fun step, args:"; val (((ip as (_,p_)):pos'), ((ptp as (pt,p), tacis):calcstate)) = 
    42                            (p, ((pt, e_pos'),[]));
    43 val pIopt = get_pblID (pt,ip);
    44 ip = ([],Res); "false";
    45 tacis; " = []";
    46 pIopt; (* = SOME ["Inverse", "Z_Transform", "SignalProcessing"]*)
    47 member op = [Pbl,Met] p_ andalso is_none (get_obj g_env pt (fst p)); "false";
    48 (*nxt_solve_ (pt,ip); "WAS isalist2list applied to NON-list 'no_meth'"
    49    THIS MEANS: replace no_meth by [no_meth] in Script.*)
    50 (*WAS val ("helpless",_) = step p ((pt, e_pos'),[]) *)
    51 (*WAS val (p,_,f,nxt,_,pt) = me nxt p [] pt; "Empty_Tac instead SubProblem";*)
    52 
    53 "----------- why not nxt = Model_Problem here ? ---------";
    54 "----------- why not nxt = Model_Problem here ? ---------";
    55 "----------- why not nxt = Model_Problem here ? ---------";
    56 val (p,_,f,nxt,_,pt) = me nxt''' p''' [] pt'''; 
    57 "~~~~~ fun me, args:"; val (((_,tac):tac'_), (p:pos'), _, (pt:ptree)) = (nxt, p, c, pt);
    58 val ("ok", (_, _, ptp)) = locatetac tac (pt,p);
    59 val (pt, p) = ptp;
    60 "~~~~~ fun step, args:"; val (((ip as (_,p_)):pos'), ((ptp as (pt,p), tacis):calcstate)) =
    61                            (p, ((pt, e_pos'),[]));
    62 val pIopt = get_pblID (pt,ip);
    63 ip = ([],Res); " = false";
    64 tacis; " = []";                                         
    65 pIopt (* = SOME ["Inverse", "Z_Transform", "SignalProcessing"]*);
    66 member op = [Pbl,Met] p_ andalso is_none (get_obj g_env pt (fst p)); " = false";
    67 (*                               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ leads into
    68 nxt_solve_, which is definitely WRONG (should be nxt_specify_ for FIND_ADD).
    69 This ERROR seems to be introduced together with ctxt, concerns Apply_Method without init_form.
    70 See TODO.txt
    71 *)
    72 
    73 "----------- fun factors_from_solution ------------------";
    74 "----------- fun factors_from_solution ------------------";
    75 "----------- fun factors_from_solution ------------------";
    76 val ctxt = ProofContext.init_global @{theory Isac};
    77 val SOME t = parseNEW ctxt "factors_from_solution [(z::real) = 1 / 2, z = -1 / 4]";
    78 val SOME (_, t') = eval_factors_from_solution "" 0 t thy;
    79 if term2str t' =
    80  "factors_from_solution [z = 1 / 2, z = -1 / 4] = (z - 1 / 2) * (z - -1 / 4)"
    81 then () else error "factors_from_solution broken";
    82 
    83 "----------- Logic.unvarify_global ----------------------";
    84 "----------- Logic.unvarify_global ----------------------";
    85 "----------- Logic.unvarify_global ----------------------";
    86 val A_var = parseNEW ctxt "A::real" |> the |> Logic.varify_global
    87 val B_var = parseNEW ctxt "B::real" |> the |> Logic.varify_global
    88 
    89 val denom1 = parseNEW ctxt "(z + -1 * (1 / 2))::real" |> the;
    90 val denom2 = parseNEW ctxt "(z + -1 * (-1 / 4))::real" |> the;
    91 
    92 val test = HOLogic.mk_binop "Groups.plus_class.plus"
    93   (HOLogic.mk_binop "Rings.inverse_class.divide" (A_var, denom1),
    94    HOLogic.mk_binop "Rings.inverse_class.divide" (B_var, denom2));
    95 
    96 if term2str test = "?A / (z + -1 * (1 / 2)) + ?B / (z + -1 * (-1 / 4))" then ()
    97   else error "HOLogic.mk_binop broken ?";
    98 
    99 (* Logic.unvarify_global test
   100 ---> exception TERM raised (line 539 of "logic.ML"): Illegal fixed variable: "z"
   101 thus we need another fun var2free in termC.sml*)
   102 
   103 "----------- eval_drop_questionmarks --------------------";
   104 "----------- eval_drop_questionmarks --------------------";
   105 "----------- eval_drop_questionmarks --------------------";
   106 if term2str test = "?A / (z + -1 * (1 / 2)) + ?B / (z + -1 * (-1 / 4))"
   107 then () else error "test setup (ctxt!) probably changed";
   108 
   109 val t = Const ("Partial_Fractions.drop_questionmarks", HOLogic.realT --> HOLogic.realT) $ test;
   110 
   111 val SOME (_, t') = eval_drop_questionmarks "drop_questionmarks" 0 t @{theory Isac};
   112 if term2str t' = "drop_questionmarks (?A / (z + -1 * (1 / 2)) + ?B / (z + -1 * (-1 / 4))) =" ^
   113   "\nA / (z + -1 * (1 / 2)) + B / (z + -1 * (-1 / 4))"
   114 then () else error "eval_drop_questionmarks broken";
   115 
   116 "----------- = me for met_partial_fraction --------------";
   117 "----------- = me for met_partial_fraction --------------";
   118 "----------- = me for met_partial_fraction --------------";
   119   val fmz =                                             
   120     ["functionTerm (3 / (z * (z - 1 / 4 + -1 / 8 * (1 / z))))", 
   121       "solveFor z", "decomposedFunction p_p"];
   122   val (dI',pI',mI') =
   123     ("Partial_Fractions", 
   124       ["partial_fraction", "rational", "simplification"],
   125       ["simplification","of_rationals","to_partial_fraction"]);
   126   val (p,_,f,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))];
   127   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   128   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   129   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   130   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   131   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   132   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   133   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   134   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   135   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   136   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   137   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   138   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   139   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   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;
   144   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   145   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   146   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   147   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   148   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   149   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   150   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   151   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   152   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   153   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   154   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   155   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   156   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   157   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   158   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   159   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   160   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   161   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   162   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   163   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   164   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   165   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   166   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   167   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   168   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   169   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   170   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   171   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   172   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   173   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   174   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   175   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   176   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   177   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   178   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   179   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   180   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   181   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   182   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   183   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   184   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   185   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   186   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   187   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   188   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   189   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   190   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   191   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   192   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   193   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   194   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   195   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   196   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   197   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   198   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   199   val (p,_,f,nxt,_,pt) = me nxt p [] pt; 
   200   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   201   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   202   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   203   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   204   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   205   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   206   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   207   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   208   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   209   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   210   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   211   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   212   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   213   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   214   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   215   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   216   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   217   val (p,_,f,nxt,_,pt) = me nxt p [] pt;
   218   val (p,_,f,nxt,_,pt) = me nxt p [] pt; 
   219 
   220 show_pt pt;
   221 if f2str f = "4 / (z - 1 / 2) + -4 / (z - -1 / 4)" then() 
   222 else error "= me .. met_partial_fraction broken";
   223 
   224 "----------- progr.vers.2: check erls for multiply_ansatz";
   225 "----------- progr.vers.2: check erls for multiply_ansatz";
   226 "----------- progr.vers.2: check erls for multiply_ansatz";
   227 (*test for outcommented 3 lines in script: is norm_Rational strong enough?*)
   228 val t = str2term "(3 / ((-1 + -2 * z + 8 * z ^^^ 2) *3/24)) = (3 / ((z - 1 / 2) * (z - -1 / 4)))";
   229 val SOME (t', _) = rewrite_set_ @{theory Isac} true ansatz_rls t;
   230 term2str t' = "3 / ((-1 + -2 * z + 8 * z ^^^ 2) * 3 / 24) =\n?A / (z - 1 / 2) + ?B / (z - -1 / 4)";
   231 
   232 val SOME (t'', _) = rewrite_set_ @{theory Isac} true multiply_ansatz t'; (*true*)
   233 term2str t'' = "(z - 1 / 2) * (z - -1 / 4) * 3 / ((-1 + -2 * z + 8 * z ^^^ 2) * 3 / 24) =\n" ^
   234   "?A * (z - -1 / 4) + ?B * (z - 1 / 2)"; (*true*)
   235 
   236 val SOME (t''', _) = rewrite_set_ @{theory Isac} true norm_Rational t'';
   237 if term2str t''' = "3 = ?A * (1 + 4 * z) / 4 + ?B * (-1 + 2 * z) / 2" then ()
   238 else error "ansatz_rls - multiply_ansatz - norm_Rational broken";
   239 
   240 (*test for outcommented 3 lines in script: empower erls for x = a*b ==> ...*)
   241 val xxx = append_rls "multiply_ansatz_erls" norm_Rational 
   242   [Calc ("HOL.eq",eval_equal "#equal_")];
   243 
   244 val multiply_ansatz = prep_rls(
   245   Rls {id = "multiply_ansatz", preconds = [], rew_ord = ("dummy_ord",dummy_ord), 
   246 	  erls = xxx,
   247 	  srls = Erls, calc = [],
   248 	  rules = 
   249 	   [Thm ("multiply_2nd_order",num_str @{thm multiply_2nd_order})
   250 	   ], 
   251 	 scr = EmptyScr}:rls);
   252 
   253 rewrite_set_ thy true xxx @{term "a+b = a+(b::real)"}; (*SOME ok*)
   254 rewrite_set_ thy true multiply_ansatz @{term "a+b = a+(b::real)"}; (*why NONE?: GOON*)
   255 
   256 "----------- progr.vers.2: improve program --------------";
   257 "----------- progr.vers.2: improve program --------------";
   258 "----------- progr.vers.2: improve program --------------";
   259 (*WN120318 stopped due to much effort with the test above*)
   260      "Script PartFracScript (f_f::real) (zzz::real) =                    " ^(*f_f: 3 / (z * (z - 1 / 4 + -1 / 8 * (1 / z)), zzz: z*)
   261      " (let f_f = Take f_f;                                       " ^
   262      "   (num_orig::real) = get_numerator f_f;                    " ^(*num_orig: 3*)
   263      "   f_f = (Rewrite_Set norm_Rational False) f_f;             " ^(*f_f: 24 / (-1 + -2 * z + 8 * z ^^^ 2)*)
   264      "   (numer::real) = get_numerator f_f;                       " ^(*numer: 24*)
   265      "   (denom::real) = get_denominator f_f;                     " ^(*denom: -1 + -2 * z + 8 * z ^^^ 2*)
   266      "   (equ::bool) = (denom = (0::real));                       " ^(*equ: -1 + -2 * z + 8 * z ^^^ 2 = 0*)
   267      "   (L_L::bool list) = (SubProblem (PolyEq',                 " ^
   268      "     [abcFormula, degree_2, polynomial, univariate, equation], " ^
   269      "     [no_met]) [BOOL equ, REAL zzz]);                       " ^(*L_L: [z = 1 / 2, z = -1 / 4]*)
   270      "   (facs::real) = factors_from_solution L_L;                " ^(*facs: (z - 1 / 2) * (z - -1 / 4)*)
   271      "   (eql::real) = Take (num_orig / facs);                    " ^(*eql: 3 / ((z - 1 / 2) * (z - -1 / 4)) *) 
   272      "   (eqr::real) = (Try (Rewrite_Set ansatz_rls False)) eql;  " ^(*eqr: ?A / (z - 1 / 2) + ?B / (z - -1 / 4)*)
   273      "   (eq::bool) = Take (eql = eqr);                           " ^(*eq:  3 / ((z - 1 / 2) * (z - -1 / 4)) = ?A / (z - 1 / 2) + ?B / (z - -1 / 4)*)
   274 (*this has been tested below by rewrite_set_
   275      "   (eeeee::bool) = Take ((num_orig / denom * denom / numer) = (num_orig / facs));" ^(**)
   276      "   (eeeee::bool) = (Rewrite_Set ansatz_rls False) eeeee;" ^(**)
   277      "   eq = (Try (Rewrite_Set multiply_ansatz False)) eeeee;         " ^(*eq: 3 = ?A * (z - -1 / 4) + ?B * (z - 1 / 2)*) 
   278 NEXT try to outcomment the very next line..*)
   279      "   eq = (Try (Rewrite_Set equival_trans False)) eeeee;         " ^(*eq: 3 = ?A * (z - -1 / 4) + ?B * (z - 1 / 2)*) 
   280      "   eq = drop_questionmarks eq;                              " ^(*eq: 3 = A * (z - -1 / 4) + B * (z - 1 / 2)*)
   281      "   (z1::real) = (rhs (NTH 1 L_L));                          " ^(*z1: 1 / 2*)
   282      "   (z2::real) = (rhs (NTH 2 L_L));                          " ^(*z2: -1 / 4*)
   283      "   (eq_a::bool) = Take eq;                                  " ^(*eq_a: 3 = A * (z - -1 / 4) + B * (z - 1 / 2)*)
   284      "   eq_a = (Substitute [zzz = z1]) eq;                       " ^(*eq_a: 3 = A * (1 / 2 - -1 / 4) + B * (1 / 2 - 1 / 2)*)
   285      "   eq_a = (Rewrite_Set norm_Rational False) eq_a;           " ^(*eq_a: 3 = 3 * A / 4*)
   286      "   (sol_a::bool list) =                                     " ^
   287      "     (SubProblem (Isac', [univariate,equation], [no_met])   " ^
   288      "     [BOOL eq_a, REAL (A::real)]);                          " ^(*sol_a: [A = 4]*)
   289      "   (a::real) = (rhs (NTH 1 sol_a));                         " ^(*a: 4*)
   290      "   (eq_b::bool) = Take eq;                                  " ^(*eq_b: 3 = A * (z - -1 / 4) + B * (z - 1 / 2)*)
   291      "   eq_b = (Substitute [zzz = z2]) eq_b;                     " ^(*eq_b: *)
   292      "   eq_b = (Rewrite_Set norm_Rational False) eq_b;           " ^(*eq_b: *)
   293      "   (sol_b::bool list) =                                     " ^
   294      "     (SubProblem (Isac', [univariate,equation], [no_met])   " ^
   295      "     [BOOL eq_b, REAL (B::real)]);                          " ^(*sol_b: [B = -4]*)
   296      "   (b::real) = (rhs (NTH 1 sol_b));                         " ^(*b: -4*)
   297      "   eqr = drop_questionmarks eqr;                            " ^(*eqr: A / (z - 1 / 2) + B / (z - -1 / 4)*)
   298      "   (pbz::real) = Take eqr;                                  " ^(*pbz: A / (z - 1 / 2) + B / (z - -1 / 4)*)
   299      "   pbz = ((Substitute [A = a, B = b]) pbz)                  " ^(*pbz: 4 / (z - 1 / 2) + -4 / (z - -1 / 4)*)
   300      " in pbz)"
   301