OKlibrary  0.2.1.6
Basic.mac
Go to the documentation of this file.
00001 /* Oliver Kullmann, 29.7.2007 (Swansea) */
00002 /* Copyright 2007, 2008, 2009 Oliver Kullmann
00003 This file is part of the OKlibrary. OKlibrary is free software; you can redistribute
00004 it and/or modify it under the terms of the GNU General Public License as published by
00005 the Free Software Foundation and included in this library; either version 3 of the
00006 License, or any later version. */
00007 
00022 oklib_include("OKlib/ComputerAlgebra/DataStructures/Lisp/Lists.mac")$
00023 
00024 
00025 /*
00026  A "branching tuple" is a list of positive numbers, which in case of length
00027  one can also include 0 (also the empty list is possible).
00028 */
00029 
00030 /* ************************
00031    * Elementary functions *
00032    ************************
00033 */
00034 
00035 bt_min(t) := lmin(t)$
00036 bt_max(t) := lmax(t)$
00037 bt_width(t) := length(t)$
00038 bt_sum(t) := apply("+", t)$
00039 
00040 /* For the concatenation of branching tuples, "append" is to be used. */
00041 
00042 /* Composition of branching tuples at position 1: */
00043 bt_composition(a,b) := if emptyp(a) then [] else append(a[1] + b, rest(a))$
00044 /* Composition of branching tuples at position i: */
00045 bt_composition_i(a,b,i) := if emptyp(a) then [] 
00046   else append(rest(a, -(length(a) - i + 1)), a[i] + b, rest(a,i))$
00047 
00048 /*
00049  A "permutation" in S_n is a list P of length n of all the elements in 
00050  {1,...,n}; it stands for the bijection i -> P[i].
00051  See ComputerAlgebra/Algebra/Lisp/Groupoids/Groups/SymmetricGroups.mac.
00052 */
00053 
00054 /* Apply a permutation (or, more generally, an index selection): */
00055 bt_apply_permutation(perm, t) := create_list(t[i], i, perm)$
00056 
00057 /* The branching-tuple-hull of a branching tuple t w.r.t. upper-bound
00058    s on the maximum, together with the set of saturated elements: */
00059 bth(t, s) := if s < bt_max(t) then [{},{}] else 
00060  block([hull : {t}, new_size : 1, old_size : 0, saturated : {}],
00061   while new_size > old_size do (
00062     old_size : new_size,
00063     for a in hull unless new_size > old_size do (
00064       hull : union(hull, permutations(a), 
00065                map(lambda([b],bt_composition(a,b)),
00066                    subset(hull,lambda([b], is(a[1]+bt_max(b) <= s))))),
00067       new_size : length(hull))),
00068   return([hull, block([mt:bt_max(t)], subset(hull, lambda([a],is(bt_min(a)+mt>s))))]))$
00069 
00070 /* Compute greedily one saturated tuple in the branching-tuple-hull: */
00071 bth_saturated_greedy(t,s) := block([m : bt_max(t), ts : sort(t)],
00072   while ts[1] + m <= s do (
00073     ts : bt_composition(ts,t),
00074     ts : sort(ts)
00075   ),
00076   return(ts))$
00077 
00078 /* The integer sequence corresponding to an integer branching tuple t,
00079    that is, the first n elements.
00080    Prerequisites: t is not empty, and also t <> [0].
00081  */
00082 int_seq_bt(t) := block([m : bt_max(t), l : length(t)], 
00083   return(buildq([m,l,t],
00084    lambda([n], block([S : create_list(1,i,1,min(n,m))],
00085      for i : m thru n-1 do
00086        S : endcons(sum(S[i-t[j]+1],j,1,l),S),
00087      return(S)))
00088    )))$
00089 /* The sequence is f(0),f(1),f(2),...,f(n-1), where
00090    f(n) := f(n-first(t)) + ... + f(m-last(t))
00091    if this expression is defined (i.e., iff n >= max(t)),
00092    while otherwise f(n) := 1.
00093    For example with fibi : int_seq_bt([1,2]) we get
00094    fibi(7) = [1, 1, 2, 3, 5, 8, 13]
00095    (note here that fib is the predefined Fibonacci-function).
00096 */
00097 
00098 
00099 /* ***************
00100    * Power means *
00101    ***************
00102 */
00103 
00104 /* n-ary arithmetic mean (while the arithmetic mean for branching tuples t
00105    is just mean(t)): */
00106 meann([t]) := mean(t)$
00107 
00108 /* Geometric mean: */
00109 oklib_plain_load(descriptive)$
00110 gmean(t) := geometric_mean(t)$
00111 /* n-ary geometric mean: */
00112 gmeann([t]) := gmean(t)$
00113 
00114 /* Harmonic mean: */
00115 hmean(t) := harmonic_mean(t)$
00116 /* n-ary harmonic mean: */
00117 hmeann([t]) := hmean(t)$
00118 
00119 /* Power means with index -inf <= p <= inf: */
00120 genmean(t, p) := if p = minf then lmin(t) elseif p = 0 then gmean(t) elseif p = inf then lmax(t) else mean(t^p)^(1/p)$
00121 /* (n+1)-ary power means: */
00122 genmeann([u]) := genmean(rest(u,-1), last(u))$
00123 /* Binary power mean: */
00124 genmean2(x,y,p) := genmeann(x,y,p)$
00125 /* Ternary power mean: */
00126 genmean3(x,y,z,p) := genmeann(x,y,z,p)$
00127 
00128 
00129 /* ***********************************************
00130    * Lower and upper bounds for the tau-function *
00131    ***********************************************
00132 */
00133 
00134 load (descriptive)$
00135 
00136 /* The lower bound given by the arithmetic mean: */
00137 tau_lo(t) := length(t)^(1/mean(t))$
00138 /* n-ary version: */
00139 tau_lon([t]) := tau_lo(t)$
00140 
00141 /* The upper bound given by the power means (defined below): */
00142 tau_up(t) := length(t)^(1/genmean(t,2-length(t)))$
00143 /* n-ary version: */
00144 tau_upn([t]) := tau_up(t)$
00145 
00146 
00147 /* ****************************
00148    * Comparison of tau-values *
00149    ****************************
00150 */
00151 
00152 /* bt_comparison(t,alpha) is pos, zero, neg if tau(t) > alpha resp. 
00153    = alpha resp. < alpha 
00154    (this decision uses only symbolic means, no numeric ones): */
00155 bt_comparison(t,alpha) := sign(chi(t,alpha) - 1)$
00156 
00157 
00158 /* *******************************************************
00159    *   Computing the tau-function                        *
00160    *******************************************************
00161 */
00162 
00163 /* The chi-functional: */
00164 chi(t,x) := bt_sum(x^(-t))$
00165 
00166 /* For fixed t, compute the Newton-step for chi(t,x)-1 = 0: */
00167 chi_newtonstep(t) := buildq([t],
00168  lambda([x], block([p : x^(-t)],
00169   return((apply("+",p)-1) / apply("+",t*p/x)))))$
00170 
00171 /* Computing the tau-function by the Newton-method: */
00172 tau_eps(t,eps) := block([l:length(t)],
00173  if l = 0 then inf elseif l = 1 then 1 else
00174  block([dxf : chi_newtonstep(t), x : tau_lo(t), dx],
00175   do (
00176     dx : dxf(x),
00177     x : x + dx,
00178     if dx <= eps then return(true)
00179   ),
00180   return(x)))$
00181 
00182 /* Given a branching tuple, compute its tau-value (supplying a standard eps):
00183 */
00184 tau(t) := tau_eps(float(t), 10^-15)$
00185 /* floatnump(tau(t)) is true */
00186 /* n-ary version (using the standard eps): */
00187 taun([v]) := tau(v)$
00188 
00189 /* Binary and ternary versions: */
00190 tau2(x,y) := taun(x,y)$
00191 tau3(x,y,z) := taun(x,y,z)$
00192 
00193 
00194 /* ********************
00195    * Higher precision *
00196    ********************
00197 */
00198 
00199 /* General remarks on floating-point representations:
00200    Use float(x) to compute a floating-point representation of x. 
00201    To use higher precision, set fpprec to the required number of
00202    decimal digits, and convert via bfloat. 
00203    Or use the following conversion function.
00204 */
00205 /* Convert a number to a number with a given precision: */
00206 hp_float(x, decimal_digits) := block([fpprec : fpprec], fpprec : decimal_digits, bfloat(x))$
00207 
00208 /* The tau-value with a given precision: */
00209 tau_hp(t, decimal_digits) := block(
00210  [fpprec : fpprec], 
00211   fpprec : decimal_digits,
00212   tau_eps(bfloat(t), 10^(-decimal_digits+1)))$
00213 
00214 
00215 /* ************************
00216    * Symbolic computation *
00217    ************************
00218 */
00219 
00220 oklib_plain_include(sqdnst)$
00221 
00222 /* Computes a symbolical solution if possible (otherwise returns false): */
00223 tau_symbolical(t) := block([l : length(t), x, sol, pos_sol, n_pos], 
00224  if l = 0 then inf elseif l = 1 then 1
00225  else (sol : map(rhs, solve(chi(t,x)-1, x)), 
00226        pos_sol : sublist_indices(sol, lambda([x], is(equal(imagpart(x),0)) and is(sign(x) = pos))),
00227        n_pos : length(pos_sol),
00228        if n_pos = 1 then 
00229          return(sqrtdenest(sqrtdenest(sol[pos_sol[1]])))
00230        elseif n_pos >= 2 then print("INCONSISTENCY with tau_symbolical for ", t)
00231        else return(false)))$
00232 /* Remark: Apply "optimize" to the result to obtain a representation with
00233    shared subterms. */
00234 
00235 /* ***********************************
00236    * Derivatives of the tau-function *
00237    ***********************************
00238 */
00239 
00240 /* The total differential for the tau-function: */
00241 Dtau(t) := block([tv : tau(t)], block([pv : tv^(-t)],
00242   return((-tv * log(tv) / apply("+", t * pv)) * pv)))$
00243 
00244 Dtaun([t]) := Dtau(t)$
00245 Dtau2(x,y) := Dtaun(x,y)$
00246 Dtau3(x,y,z) := Dtaun(x,y,z)$
00247 
00248 /* Symbolic differentiation of tau2: */
00249 Dtau2s(x,y) := block( [ tv : tau2s(x,y)], block( [pv : tv^(-[x,y]) ] ,
00250   return((- tv * log(tv) / apply("+", [x,y] * pv)) * pv)))$
00251 
00252 gradef(tau2s(x,y), Dtau2s(x,y)[1], Dtau2s(x,y)[2])$
00253 
00254 
00255 /* ****************************************
00256    * The induced probability distribution *
00257    ****************************************
00258 */
00259 
00260 /* The probability distribution derived from a non-empty branching tuple: */
00261 tauprob(t) := if length(t)=1 then [1] else
00262  block([ts : t / lmin(t)], tau(ts)^(-ts))$
00263 tauprobn([t]) := tauprob(t)$
00264 tauprob2(x,y) := tauprobn(x,y)$
00265 tauprob3(x,y,z) := tauprobn(x,y,z)$
00266 
00267 /* Higher precision: */
00268 tauprob_hp(t, decimal_digits) := if length(t)=1 then [1] else
00269  block([fpprec : fpprec], 
00270   fpprec : decimal_digits, t : bfloat(t),
00271   block([ts : t / lmin(t)], tau_hp(ts,decimal_digits)^(-ts)))$
00272 
00273 /* Symbolical computation (if possible): */
00274 tauprob_symbolical(t) := tau_symbolical(t)^(-t)$
00275 /* Applying factor(expand(r)), radcan(r) or fullratsimp(r) to the result r 
00276    might help simplifying. */
00277 
00278 /* The canonical inversion for tuples p with bt_sum(p) = 1: */
00279 invtauprob(p) := -map(log,p)$
00280 /* tauprob(invtprop(p)) = p and invtprop(tauprob(t)) = t * log(tau(t)) */
00281 
00282 
00283 /* ****************
00284    * The tau-mean *
00285    ****************
00286 */
00287 
00288 /* Auxiliary function to test for zero: */
00289 zerop(x) := is(equal(x,0))$
00290 /* Auxiliary function to test for +infinity: */
00291 infp(x) := is(equal(x,inf))$
00292 
00293 /* The basic tau-mean: */
00294 taumeanb(t) := float(log(length(t))) / log(tau(t))$
00295 taumeanb2(x,y) := taumeanb([x,y]);
00296 taumeanb3(x,y,z) := taumeanb([x,y,z]);
00297 /* The tau-mean for tuples allowing zero-values: */
00298 taumean(t) := if length(t) <= 1 then inf 
00299   elseif some_s(zerop, t) then 0 
00300   else float(log(length(t))) / log(tau(t))$
00301 /* The n-ary tau-mean: */
00302 taumeann([t]) := taumean(t)$
00303 /* The binary tau-mean: */
00304 taumean2(x,y) := taumeann(x,y)$
00305 /* The ternary tau-mean: */
00306 taumean3(x,y,z) := taumeann(x,y,z)$
00307 
00308 /* Tau-mean additionally also allowing inf-values: */
00309 taumean_inf(t) := if some_s(infp, t) then 
00310     block([t : delete(inf, t)], if t = [0] then unknown else taumean(t))
00311   else taumean(t)$
00312 taumean_infn([t]) := taumean_inf(t)$
00313 
00314 
00315 /* *********************************************
00316    *   Convexity considerations: line versions *
00317    *********************************************
00318 */
00319 
00320 
00321 /* Considerung for a branching tuple t0 and a direction d0 (a real vector of
00322    the same length as t) the function x -> tau(t0 + x * d0): */
00323 
00324 tau_line(x) := (local(x), tau(t0 + x * d0))$
00325 logtau_line(x) := (local(x), log(tau(t0 + x * d0)))$
00326 loglogtau_line(x) := (local(x), log(logtau_line(x)))$
00327 tauprob_line(x) := (local(x), tauprob(t0 + x * d0)[i0])$
00328 taumean_line(x) := (local(x), taumean(t0 + x * d0))$
00329 /* assume(length(t0) >= 1), assume(t0[i] > 0) */
00330 taumeandiff_line(x) := (local(x, t), 
00331   block([t : t0 + x * d0], taumean(t) - genmean(t, p0))
00332 )$
00333 
00334 /* Returning the interval bounds (for the open interval): */
00335 tau_line_dom() := block(
00336   [pos : sublist_indices(d0, lambda([x], x > 0)), neg : sublist_indices(d0, lambda([x], x < 0))],
00337   [lmax(- part(t0, pos) / part(d0,pos)), lmin(- part(t0, neg) / part(d0, neg))]
00338 )$
00339 
00340 
00341 /* The chi-function on a line (but chi(t,x) only for x > 1): */
00342 chi_line(x) := (local(x), chi(t0 + x * d0, a0 + x * b0))$
00343 /* assume(length(t0) >= 2), assume(t0[i] > 0), assume(a0 > 1) */
00344 
00345 chi_line_dom() :=
00346  block([ts : append(t0, [a0 - 1]), ds : append(d0, [b0])], 
00347    block(
00348     [pos : sublist_indices(ds, lambda([x], x > 0)), neg : sublist_indices(ds, lambda([x], x < 0))],
00349     [lmax(- part(ts, pos) / part(ds,pos)), lmin(- part(ts, neg) / part(ds, neg))]
00350 ))$
00351 
00352 chin([u]) := chi(rest(u,-1), last(u))$
00353 chi2(t1,t2,x) := chin(t1,t2,x)$
00354 chi3(t1,t2,t3,x) := chin(t1,t2,t3,x)$
00355 
00356 hess_chi2 : hessian(chi2(t1,t2,x),[t1,t2,x])$
00357 /* if for example x has a value, then apply remvalue(x), or, better, kill(x) */
00358 
00359 eigen_chi2(a,b,c) := eigens_by_jacobi(at(hess_chi2, [t1=a,t2=b,x=c]))$
00360 /* a, b > 0, c > 1 */
00361 test_posdef_chi2(a,b,c) := lmin(eigen_chi2(a,b,c)[1])$
00362 
00363 
00364 /* ************************************
00365    * Investigations on approximations *
00366    ************************************
00367 */
00368 
00369 /* Consider the bounds
00370      genmean(t,k-2) <= taumeanb(t) <= mean(t).
00371 */
00372 
00373 /* For k = 2 the lower bound is always better as an approximation than
00374    the upper bound.
00375    Here we investigate how this can be quantified. 
00376 */
00377 /* For k >= 3 this is not always the case (example [5.0,5.0,0.835]),
00378    but it appears that still the lower bound yields a better approximation:
00379     - perhaps regarding the integral. ?
00380 */
00381 /* It is now known yet what happens for k > 3. */
00382 
00383 dltau_2(x,y) := taumeanb2(x,y) - float(gmean([x,y]))$
00384 dutau_2(x,y) := mean([x,y]) - taumeanb2(x,y)$
00385 dtau_2(x,y) := dutau_2(x,y) - dltau_2(x,y)$
00386 /* For all x,y > 0 we have dtau_2(x,y) >= 0 (equality holds
00387    iff x = y). */
00388 
00389 dltau_3(x,y,z) := taumeanb3(x,y,z) - float(hmean([x,y,z]))$
00390 dutau_3(x,y,z) := mean([x,y,z]) - taumeanb3(x,y,z)$
00391 dtau_3(x,y,z) := dutau_3(x,y,z) - dltau_3(x,y,z)$
00392 
00393