It's not clear to me what problem your proposal is intending to solve. In what way is it different from the current evalf algorithm/results?
-- Oscar On Wed, 11 Mar 2020 at 22:48, Shubham thorat <sthorat...@gmail.com> wrote: > > This is how I have divided the tasks: > The algorithm is defined in this paper: > https://herbie.uwplse.org/pldi15-paper.pdf > > Error Calculation on sampling points > > Find the best-accepted level of precision so we calculate the actual > correctly up to 64 bits by increasing the precision until we get constant up > to 64 bits. > Calculate the value using float( hardware precision). > Compare real and float value by calculating E(x,y) = log(2,z), z = number of > floating-point between real and approximate answers. > averaging these differences over the sampling to see how accurate the > expression is. > > Pick Candidate > > Pick candidates (subexpression) from the table and apply error calculation as > well as a local error at each location on the sampled points. > the database will have a set of rewrite rules like commutativity, > associativity, distributivity, (x + y) = (x**2 - y**2)/(x - y), (x - y) = > (x**3 - y**3)/(x**2 + y**2 + x*y), x = exp(log( x )) etc.. > > Recursive- rewrite > > Rewrite candidates using the database of rules and simplify to cancel terms. > Recursively repeat the algorithm on the best subexpression. > > Series Expansion > > Finding the series expansion of expressions to remove error near 0 and > infinity. > > Candidate Tree > > Only keep those candidates that give the best accuracy on at least one > location. > On every iteration of the outer loop, the algorithm chooses the program from > the table and uses it to find new candidates, every program is used once. > Candidate programs are also saved as a pair of maps that are tied with the > location that they are best at. > removing candidates if more than one candidate is giving the same results > based on their results at other locations. > Before the series approximation step, we will use the set cover approximation > algorithm to prune candidates to have a minimal set. > > Get Piecewise solutions > > Split is found using dynamic programming and later refined using binary > search. > > > These are the functions: > > > Definition-main(program) : > > points := sample-inputs(program) > > exacts := evaluate-exact(program, points) > > table := make-candidate-table(simplify(program)) > > repeat N times > > candidate:= pick-candidate(table) > > locations := sort-by-local-error(all-locations(candidate)) > > locations.take(M ) > > rewritten := recursive-rewrite(candidate, locations) > > new-candidates := simplify-each(rewritten) > > table.add(new-candidates) > > approximated := series-expansion(candidate, locations) > > table.add(approximated) > > return infer-regimes(table).as-program > > > > Definition local-error(expr, points) : > > for point ∈ points : > > args := evaluate-exact(expr.children) > > exact-ans := F(expr.operation.apply-exact(args)) > > approx-ans := expr.operation.apply-approx(F(args)) > > accumulate E(exact-ans, approx-ans) > > > > Definition recursive-rewrite(expr, target) : > > select input > > output from RULES > > where input.head = expr.operator ∧ output.head = target.head > > for (subexpr, subpattern) ∈ zip(expr.children, input.children) : > > if ¬matches(subexpr, subpattern) : > > recursive-rewrite(subexpr, subpattern) > > where matches(expr, input) > > expr.rewrite(input -> output) > > > > Definition infer-regimes(candidates, points) : > > for x i ∈ points : > > best-split 0 [x i ] = [regime(best-candidate(−∞, x i ), −∞, x i )] > > for n ∈ N until best-split n+1 = best-split n : > > for x i ∈ points ∪ {∞} : > > for x j ∈ points, x j < x i : > > extra-regime := regime(best-candidate(x j , x i ), x i , x j ) > > option[x j ] := best-split[x j ] ++ [extra-regime] > > best-split n+1 [x i ] := lowest-error(option) > > if best-split n [x i ].error − 1 ≤ best-split n+1 [x i ].error : > > best-split n+1 [x i ] := best-split n [x i ] > > split:= best-split ∗ [∞] > > split.refine-by-binary-search() > > return split > > > > I have written a basic brute force code. > from sympy import * > import numpy as np > > x = Symbol("x") > y = Symbol("y") > z = Symbol("z") > > > #points > maxi = 1000000000000000000000000000000 > mini = -1000000000000000000000000000000 > step = (maxi-mini)/256 > > start = mini > points = [] > > for i in range(0,256): > points.append(start) > start += step > > #calculate error > def calc_error(expr,point): > from mpmath import mp, mpf > mp.dps = 1000 > symb = expr.atoms(Symbol) > unq_sym = len(symb) > > > subst_sym = [] > subst_mpf = [] > > i=0 > for sym in symb: > subst_sym.append((x,point)) > #subst_sym.append((x,300)) > #subst_sym.append((z,400)) > subst_mpf.append( ( sym,mpf(point) ) ) > i = i+1 > > ans1 = expr.subs(subst_sym) > ans2 = expr.subs(subst_mpf) > return ans1,ans2 > > #replacement functions > #database > rule = [] > rule.append(lambda exp: (exp.args[0]**3 + exp.args[1]**3)/(exp.args[0]**2 + > exp.args[1]**2 - exp.args[1]*exp.args[0]) if type(exp) ==Add and > len(exp.args)==2 and (type(exp.args[0]) != Symbol and type(exp.args[1]) != > Symbol) else False) > rule.append(lambda exp: (exp.args[0]**2 - > exp.args[1]**2)/(exp.args[0]-exp.args[1]) if type(exp) ==Add and > len(exp.args)==2 and (type(exp.args[0]) != Symbol and type(exp.args[1]) != > Symbol)else False) > #rule.append(lambda expr: exp(log(expr))) > > #rule1 = lambdify([x, y, (x**2 - y**2)/(x+y)) > > expr = (x+1)**0.5 - x**0.5 > #expr = x+1-x #+ (1/(x-1)) - (2/x) > #expr = (-y + (-4*x*z + y**2)**0.5)/(2*x) > #expr = (x+1)**(1/3) - x**(1/3) > #expr = (x+1)**0.25 - x**0.25 > expr = simplify(expr) > main_expr = expr > temp_expr = expr > pprint(rule[0](expr)) > print(rule[1](expr)) > > all_func = [] > mapper = dict() > > def pre(expr): > for fi in rule: > print("see ",expr) > for pnt in points: > for fi in rule: > k = fi(expr) > if k!=False: > ans1, ans2 = calc_error(expr, pnt) > update_func = main_expr.subs(expr, k) > ans3, ans4 = calc_error(update_func, pnt) > diff1 = ans2 - ans1 > diff2 = ans4 - ans3 > #print("istrue ", ans4==ans2) > #print("check1 ", ans4) > #print("check2 ", ans2) > #print("difference ", diff1, " ", diff2) > try: > if abs(diff2) <= abs(diff1): #(diff2) <= (diff1) > all_func.append(update_func) > print("inn ", pnt) > try: > lister = mapper[update_func] > lister.append(pnt) > mapper[update_func] = lister > except: > mapper[update_func] = [pnt] > except: > print("failed ",pnt) > #print(abs(diff1), " ", abs(diff2)) > #gprint(isinstance(diff2, complex)) > #print(update_func) > #exit() > #all_func.append(simplify(main_expr.subs(expr,k))) > print(expr, " space ",) > for arg in expr.args: > pre(arg) > > if len(all_func)==0: > all_func.append(simplify(main_expr)) > > pre(expr) > #print(srepr(expr)) > #pprint(set(all_func)) > #print("mapper") > #print(mapper) > values = list(set(all_func)) > i=0 > for expr in values: > values[i] = expr.replace( lambda expr: type(expr)== Pow and expr.args[1]==1, > lambda expr: Mul(expr.args[0],1)) > i = i+1 > pprint(values)Enter code here... > > Output: > > 1.5 1.5 > - x + (x + 1) 1 > 0.5 0.5 > ────────────── , ────────── , - x + (x + 1) > 0.5 0.5 0.5 0.5 > x ⋅ (x + 1) + 2⋅x + 1 x + (x + 1) > > > > > please tell me your thoughts on this, any addition, deletion of algorithmic > improvement. > > -- > You received this message because you are subscribed to the Google Groups > "sympy" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to sympy+unsubscr...@googlegroups.com. > To view this discussion on the web visit > https://groups.google.com/d/msgid/sympy/f94a7d30-dc4c-4ba0-a333-e67ff1b63513%40googlegroups.com. -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to sympy+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxRa3LRO1iVCs3uL5gCKV-h9ciYUDz4XpKOURJwkmohOmw%40mail.gmail.com.
