text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: work_mgga_c *)

$define gga_x_b88_params
$include "gga_x_b88.mpl"

cab := 0.63:
css := 0.96:

b88_cab := (rs, z, xs0, xs1) ->
  - 0.8 * (1 - z^2)/4 * n_total(rs)
  * b88_zab(cab, f_b88, rs, z, xs0, xs1)^2 * (
    1 - log(1 + b88_zab(cab, f_b88, rs, z, xs0, xs1))
      / b88_zab(cab, f_b88, rs, z, xs0, xs1)
    ):

b88_css := (rs, z, xs, ts) ->
  - 0.01 * (1 + z)/2 * n_spin(rs, z)^(5/3) * 2*ts * Fermi_D(xs, ts)
  * b88_zss(css, f_b88, rs, z, xs)^4 * (
    1 - 2*log(1 + b88_zss(css, f_b88, rs, z, xs)/2)
      / b88_zss(css, f_b88, rs, z, xs)
    ):

f_b88_c := (rs, z, xs0, xs1, ts0, ts1) ->
  + b88_cab(rs,  z, xs0, xs1)
  + b88_css(rs,  z, xs0, ts0)
  + b88_css(rs, -z, xs1, ts1):

f := (rs, z, xt, xs0, xs1, ts0, ts1, us0, us1) ->
  f_b88_c(rs, z, xs0, xs1, ts0, ts1):


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

op_a1 := 1.5214:
op_a2 := 0.5764:
op_b1 := 1.1284:
op_b2 := 0.3183:

op_beta := (rs, z, xs0, xs1) ->
  op_qab/b88_zab(1, op_f, rs, z, xs0, xs1):

f_op := (rs, z, xt, xs0, xs1) ->
  - (1 - z^2)*n_total(rs)/4.0
  * (op_a1*op_beta(rs, z, xs0, xs1) + op_a2)
  / (op_beta(rs, z, xs0, xs1)^4 + op_b1*op_beta(rs, z, xs0, xs1)^3 + op_b2*op_beta(rs, z, xs0, xs1)^2):

f := (rs, z, xt, xs0, xs1) ->
  f_op(rs, z, xt, xs0, xs1):

ext: mpl

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)

params_a_beta  := 0:
params_a_alpha := 0.8:
$include "gga_c_zpbeint.mpl"

(* parameters for beta of SG4, b0 = 3 mu^MGE2/pi^2, mu^MGE2 = 0.26 *)
sg4_b0    := 3*0.262/Pi^2:
sg4_sigma := 0.07:

(* we redefine beta here *)
mbeta := (rs, t) -> sg4_b0 + sg4_sigma*t*(1 - exp(-rs^2)):


text:
# real 3 metric for BSSN in spherical symmetry.
Ndim_ := 3:
x1_ := x:
x2_ := theta:
x3_ := phi:
complex_ := {}:
g11_ := A(x)*exp(4*phi(x)):
g12_ := 0:
g13_ := 0:
g22_ := x^2*B(x)*exp(4*phi(x)):
g23_ := 0:
g33_ := x^2*B(x)*(sin(theta))^2*exp(4*phi(x)):


text:
with(DifferentialAlgebra):
ring_diff := DifferentialRing(blocks = [[x3, x1, x2, x0], [y1] ], derivations = [t]):
sys := [
diff(x3(t), t) - (-b2*x3(t) + b2*x2(t)),
diff(x1(t), t) - (a2*x0(t) - a2*x1(t)),
diff(x2(t), t) - ((ka*kc*b1*x3(t) - ka*kc*b1*x2(t) + ka*b1*x0(t)*x3(t) - ka*b1*x0(t)*x2(t) - n*kc*x2(t) + kc*b1*x3(t)*x2(t) - kc*b1*x2(t)^2) / (ka*kc + ka*x0(t) + kc*x2(t))),
diff(x0(t), t) - ((-ka*n*x0(t) - ka*kc*a1*x0(t) + ka*kc*a1*x1(t) - ka*a1*x0(t)^2 + ka*a1*x0(t)*x1(t) - kc*a1*x0(t)*x2(t) + kc*a1*x1(t)*x2(t)) / (ka*kc + ka*x0(t) + kc*x2(t))),
y1(t) - (x0(t))
];
res := CodeTools[CPUTime](RosenfeldGroebner(sys, ring_diff, singsol=none));


text:
use DocumentTools in

SetProperty("sigma", value, ""):
exname := GetProperty("example_box", value):
# cleanup
SetProperty("LogAreaSIAN", value, ""):
SetProperty("LogAreaSE", value, ""):
SetProperty("LogAreaME", value, ""):
SetProperty("GlobalParams1", expression, NULL):
SetProperty("LocalParams1", expression, NULL):
SetProperty("NoIDParams1", expression, NULL):
SetProperty("being_refined", caption, ""):
SetProperty("Bound", expression, NULL):
SetProperty("MultiFunctions", expression, NULL):
SetProperty("SingleFunctions", expression, NULL):
SetProperty("RunningTimeSIAN", value, ""):
SetProperty("RunningTimeMulti", value, ""):
SetProperty("RunningTimeSingle", value, ""):
SetProperty("Parameters", value, ""):

SetProperty("sigma", value, exname):
if exname = "Custom" then
	SetProperty("sigma", value, ""):
else
	SetProperty("sigma", value,cat(op(examples[exname][2])),'refresh'):
fi:
end use;

max_stars_repo_name: iliailmer/sian-web-app

max_issues_repo_name: iliailmer/sian-web-app

max_forks_repo_name: iliailmer/sian-web-app


text:
# Date: Thu, 24 Mar 2016 18:55:55 -0400
# From: Chung-chieh Shan <ccshan@indiana.edu>
# To: Jacques Carette <carette@mcmaster.ca>
# Subject: Re: partial implementation
#
# Great, let's say a KB means a set of tuples -- in other words, a n-ary
# relation.  (I'm going to treat all variables as parameters for now.)
# So "empty" means the singleton set containing unit, and "gen*" means
# to increment the dimensionality n by 1, such as taking the Cartesian
# product of a relation with an interval.  And "assert" makes a subset.
#
# [...]  let me define subtraction semantically.  Subtraction means
# coming up with a (hopefully compact) sequence of "gen*" and "assert"
# operations that take you from one KB to another.  Of course, for that
# to be possible, the source KB must be a superset of a projection of the
# destination KB.  (Here by projection I mean in the relational database
# sense.)
# TODO: decide if this should be kept

# KB, short for knowledge base, is an 'abstract' representation
# of constraints on bound variables. Abstractly, it is a (boolean) conjunction
# of clauses, each of which is of one of the following forms:
#  - Introduce(n::Name,t::HakaruType)
#      introduces a new variable "n : t" which scopes over
#      the 'rest' of the KB (to the right (?))
#  - Let(n::Name,e::Hakaru)
#      introduces a new variable "n" which scopes over the rest
#      of the KB whose value is precisely "e"
#  - Constrain(c::Constraint) where
#      where Constraint is essentially a boolean-valued expression
#      on Hakaru terms, containing equalities, inequalities,
#      and boolean algebra on such
KB := module ()
  option package;
  local
     # The 'constructors' of a KB, which should not be used externally to the module
     KB, Introduce, Let, Constrain,

     # Experimental KB which allows one 'level' of a KB to be a partition along
     # a particular dimension/variable.  Under that split, what is below are other
     # layers of SplitKB.  A KB is a degenerate case of a KB with no split.
     SplitKB,

     # Some sort of hideous hack to get built-in maple
     # functions (assuming,eval,is...) to work with
     # Hakaru types
     chilled,

     ModuleLoad, ModuleUnload,

     # Various utilities
     t_intro, t_lo, t_hi, log_metric,
     boolean_if, coalesce_bounds, htype_to_property, bad_assumption, bad_assumption_pw

     ;
  export
     # Functions which build up KBs from KBs and other pieces
     #  typically ensuring that internal invariants are upheld
     # (i.e. 'smart' constructors)
      empty, genLebesgue, genType, genSummation, genIntVar, genLet,
      assert, assert_deny, assert_mb, assert_deny_mb, build_kb,

     # for debugging
     build_unsafely,

     chill, warm,

     # Negation of 'Constrain' atoms, that is, equality and
     # inequality constraints
     negated_relation, negate_rel, negate_rels,

     # "kb0 - kb1" - that is, kb0 without the knowledge of kb1
     kb_subtract,

     # kb_entails(kb,cond) = "kb => cond"
     kb_entails,

     # Simplify a Hakaru term assuming the knowledge of the kb
     # variants do different things in case of simplification error
     # (which should really only occur when the KB contains a contradicition)
     kb_assuming_mb, simplify_assuming, simplify_assuming_mb, simplify_assuming_f,

     # Gets the most refined (see refine_given) type of a given name under the
     # assumptions of the KB; & convert such a type to a range.
     getType,

     # Various 'views' of the KB, in that they take a KB and produce something
     # which is somehow 'representative' of the KB
     kb_to_variables, kb_to_assumptions, kb_to_constraints, kb_to_equations, kb_piecewise, kb_Partition,

     kb_atom_to_assumptions,

     # Various utilities ...
     list_of_mul, for_poly, range_of_HInt, eval_kb, kb_is_false, try_improve_exp,

     # Types corresponding to the constructor forms of the 'atoms' of KBs
     t_kb_Introduce, t_kb_Let, t_kb_Bound, t_kb_Constrain;
  global
     # The type of KBs. It is 'assumed' that things 'with this type'
     # are actually KBs in the proper form
     t_kb, t_kb_atom, t_kb_atoms,

     # Things which should produce a KB, but sometimes don't return this
     # expression
     NotAKB,

     # kb_LMS produces NoSol when it doesn't find a solution.
     NoSol,

     # Some silly things that KB must do to appease
     # Maple when using Maple functions to work with
     # Hakaru 'terms'.
     `expand/product`, `simplify/int/simplify`,
     `product/indef/indef`, `convert/Beta`;
  uses Hakaru, SolveTools:-Inequality ;

  # The type of introductions about which we want to pass some useful
  # information to Maple.
  t_intro := 'Introduce(name, And(t_type, Not(specfunc({HMeasure,HFunction}))))';

  # Low and high bounds (?)
  t_lo    := 'identical(`>`,`>=`)';
  t_hi    := 'identical(`<`,`<=`)';

  # The 'constructor' forms of KB
  t_kb_Introduce := 'Introduce(name, anything)';
  t_kb_Let       := 'Let(name, anything)';
  t_kb_Bound     := 'Bound(name, anything, anything)';
  t_kb_Constrain := 'Constrain(anything)';

  # The empty KB means "true".
  empty := KB();

  # The false KB
  NotAKB := proc($)
      'procname'()
  end proc;

  # Check if a kb which might be NotAKB is indeed NotAKB
  kb_is_false := proc(mbkb :: t_kb_mb, $) evalb(mbkb :: t_not_a_kb) end proc;

  # A smart constructor for introducing Lebesgue integration variables (?)
  #   genLebesgue(var,lo,hi,kb) =
  #     "KB(Introduce(x::AlmostEveryReal(x>lo,x<hi))
  #        ,kb)"
  genLebesgue := genIntVar(proc(lo,hi,$) `AlmostEveryReal`(Bound(`>`,lo), Bound(`<`, hi)) end proc);

  # Another type of integration variable
  genSummation := genIntVar(proc(lo,hi,$) `HInt`(Bound(`>=`,lo), Bound(`<=`, hi)) end proc);

  # A smart constructor for 'integration' (of which summation is a variety)
  # variables.
  genIntVar := proc (kind,$) proc(xx::name, lo, hi, kb::t_kb)
      genType(xx, kind(lo,hi), kb, _rest);
  end proc; end proc;

  # A smart constructor for type introductions. ensures name binding
  # is done correctly.
  #   genType(var,type,kb) =
  #     "KB(Introduce(x::type),kb)"
  genType := proc(xx::name, tt::t_type, kb::t_kb)
    # A variable created using genType is a parameter, in the sense that its
    # value is completely respected
    local x, t;
    x := `if`(depends([tt,kb,_rest], xx), gensym(xx), xx);
    t := subsindets(tt, identical(Bound(`>` , -infinity),
                                  Bound(`>=`, -infinity),
                                  Bound(`<` ,  infinity),
                                  Bound(`<=`,  infinity)), _->NULL);
    x, KB(Introduce(x, t), op(kb));
  end proc;

  # A smart constructor for 'let' introductions. ensures name binding
  # is done correctly.
  #   genLet(var,expr,kb) =
  #     "KB(Let(var=expr),kb)"
  genLet := proc(xx::name, e, kb::t_kb)
    # A let-binding, assigning a known value to a new variable
    local x, t;
    x := `if`(depends([e,kb,_rest], xx), gensym(xx), xx);
    x, KB(Let(x, e), op(kb));
  end proc;

  #Simplistic negation of relations. Used by Hakaru:-flatten_piecewise.
  negated_relation:= table([`<`, `<=`, `=`, `<>`] =~ [`>=`, `>`, `<>`, `=`]);

  # Takes the bool type (true/false) to mean universal and empty relations respectively.
  # i.e. negate R, where R is an 'atomic' relation of a KB
  negate_rel:= proc(R::t_kb_atom, $)::t_kb_atom;
      if R :: truefalse then
        not R
      elif R :: relation then
          negated_relation[op(0,R)](op(R));
      else
          # This is to appease 'piecewise', which won't be happy with Not
          # However, KB doesn't really care - it's already expecting {Not,not}
          # 'Technically' this is a KB 'constructor'!
          not(R);
      end if;
  end proc;

  negate_rels := proc(e, $)
      subsindets(e, { specfunc(relation, `Not`), `not`(relation) }, negate_rel@op );
  end proc;

  # Builds a kb from a list of atoms - simply foldr(assert,empty,as) except
  # and extra check is (optionally) done for the resulting KB to be valid.
  build_kb := proc(as::t_kb_atoms, shouldBeValid::{identical(false),string} := false, initKb::t_kb := empty, $)
      local
      kb := initKb;

      kb := foldr(assert,kb, op(as));
      if shouldBeValid :: string then
          ASSERT(type(kb,t_kb), sprintf("%s (in build_kb): KB contains a contradiction.", shouldBeValid));
      end if;
      kb
  end proc;

  assert_mb := proc(b::t_kb_atom, mkb::t_kb_mb, $)
      if mkb :: t_kb then
          assert(b,mkb);
      else
          NotAKB();
      end if;
  end proc;

  assert_deny_mb := proc(bb0::t_kb_atom, pol::identical(true,false), kb::t_kb_mb, $)
      if kb :: t_kb then
          assert_deny(bb0,pol,kb);
      else
          NotAKB();
      end if;
  end proc;

  # Like assert_deny, except does not accept a boolean
  # parameter to indicate negation, and evaluates
  # (using Maple's eval, anything can happen!) the
  # new conjunct under the derived knowledge of the KB
  assert := proc(b::t_kb_atom, kb::t_kb, $)
    assert_deny(foldl(eval, b, op(kb_to_equations(kb))), true, kb)
  end proc;

  # Implements the assert_deny function, which inserts either
  # a constraint or its negation into a KB, essentially, a
  # 'smart constructor'.
  # Also implements a second way to call it, 'part', which will
  # return a SplitKB instead.
  assert_deny := module ()
   export ModuleApply;
   local t_if_and_or_of, t_not, t_constraint_flipped, bound_simp, not_bound_simp, postproc_for_solve,
         refine_given, t_bound_on, simplify_in_context, expr_indp_errMsg, rel_coulditbe;

   # Either And or Or type, chosen by boolean pol
   t_if_and_or_of := proc(pol,$)
       `if`(pol, '{specfunc(anything, And), `and`}', '{specfunc(anything, Or ), `or` }')
   end proc;

   # The 'type' of `not(..)` statements
   t_not := '{specfunc(anything, Not), `not`}';

   # The type representing equalities
   # between something which is neither a name nor 'size' applied to a name
   # and another thing which is a name which is neither constant nor undefined
   t_constraint_flipped := 'Not({name, size(name)}) = Name';

   # The 'type' representing bounds on `x', where `x' is a name
   t_bound_on := proc(x,$) And(relation, Or(anyop(identical(x), freeof(x)), anyop(freeof(x), identical(x)))) end proc;

   # Given (TODO: add these types to the function(?))
   #   k  :: HakaruType
   #   kb :: t_kb
   #   x  :: name
   #   c  :: type
   # a bound name "x" of type "k", produces a 'refinement' of the
   # type "k", which is a type more specific than "k" at which
   # the name "x" is well-typed given the knowledge in the KB,
   # i.e. the conjunction of "x :: $result" with "kb" gives
   # "x :: k" inside the KB
   #
   # Note that if the list is non-empty, the refined knowledge
   #   is given by the first element of the list, and if it
   #   is empty, no refinement is possible.
   #
   # Original enlightenment:
   #   hakaru-dev/hakaru/commit/02b9335669e00921a57c3d2a65a1f5f9f6162aa4
   refine_given := proc(k,kb,x,c,$)
       [ op(map2(subsop
                , 1=NULL
                , select(type, kb, Bound(identical(x), c, anything))
                )
           )
       , op(select(type, k , Bound(              c, anything)) )
       ]
   end proc;

   # Performs simplification(?) in case something of the form `t_bound_on` is
   # found
   # This function signals it has failed to find a result with `FAIL`
   bound_simp := proc(b,x,k,kb,pol,as,$)
      local e, rel, c, kb1,ch;
      # b is a bound on x, so compare it against the current bound on x.
      # First, express `if`(pol,b,Not(b)) as rel(x,e)
      rel := op(0,b);

      # Account for the symmetric cases of `x' being on the left or right
      # hand side; the predicate is modified(?) in case of a 'flip'
      if x = lhs(b) then
        e := rhs(b);
      else #x=rhs(b)
        e   := lhs(b);
        rel := subs({`<`=`>`, `<=`=`>=`}, rel);
      end if;

      # Change relations to their negations if `pol = false'
      if not pol then
        # A 'table' giving the negations of relations
        rel := subs({`<`=`>=`, `<=`=`>`, `>`=`<=`, `>=`=`<`, `=`=`<>`, `<>`=`=`}, rel);
      end if;

      # Make inequalities loose over Real
      # Warning: other code relies on this!!!
      if k :: 'specfunc(AlmostEveryReal)' then
        # A 'table' giving the relaxation of strict relations
        rel := subs({`<=`=`<`, `>=`=`>`}, rel);
      end if;

      if rel = `=` then
        # To assert that x=e, it's not enough to supersede the Introduce
        # binding for x with a Let binding.
        kb1 := KB(Bound(x,`=`,e), op(kb));
        # We also need to assert that e is in bounds for x.
        for c in t_lo, t_hi do
          c := refine_given(k,kb,x,c);
          if nops(c)>0 then
            kb1 := assert_deny_mb(op([1,1],c)(e,op([1,2],c)), true, kb1)
          end if
        end do;
        return kb1
      end if;

      ch := chill(e);
      if rel = `<>` then
        # Refine <> to > or < if possible.
        if   is(x<=ch) assuming op(as) then rel := `<`
        elif is(x>=ch) assuming op(as) then rel := `>`
        else return KB(Constrain(x<>e), op(kb)) end if
      end if;

      # Strengthen strict inequality on integer variable.
      if op(0,k) = HInt then
        if rel = `>` then
          rel := `>=`;
          ch  := floor(ch)+1 assuming op(as);
          e   := warm(ch)
        elif rel = `<` then
          rel := `<=`;
          ch  := ceil(ch)-1 assuming op(as);
          e   := warm(ch)
        end if
      end if;

      # Look up the current bound on x, if any.
      c := `if`(rel :: t_lo, t_lo, t_hi);
      c := refine_given(k,kb,x,c);

      # chill but also unwraps `c' (?)
      if nops(c) > 0 then c := chill(op(1,c)) end if;

      # Compare the new bound rel        (x,e          )
      # against the old bound op([1,1],c)(x,op([1,2],c))
      if e = `if`(rel :: t_lo, -infinity, infinity)
          or nops(c)>0 and (is(rel(y,ch)) assuming
                              (op(1,c)(y,op(2,c)),
                               y::htype_to_property(k), op(as)) ) then
          # The old bound renders the new bound superfluous.
          return kb

      # "assuming" used to be in a try which would
      # cause the return to never be reached if it threw, but now it
      # produces FAIL instead - and `_x or FAIL = FAIL'
      elif nops(c)=0 or (is(op(1,c)(y,op(2,c))) assuming
                             (rel(y,ch),
                              y::htype_to_property(k), op(as)) ) then
          # The new bound supersedes the old bound.
          return KB(Bound(x,rel,e), op(kb))
      end if;

      FAIL; # No simplification could be done
   end proc;

   # Simplification when the `:: t_bound_on' predicate is false
   # note: k is ignored, but this makes the API the same as
   # bound_simp
   not_bound_simp := proc(b,x,k,kb,pol,as,$)
     local c, bad;
     if _Env_HakaruSolve=false then return FAIL; end if;
     if x::relation then
       userinfo(6, 'KB', printf("Chose not to solve %a\n",x));
       return FAIL;
     end if;

     # don't even try to solve bad cases, we might get a RootOf !
     bad := select(depends, indets(b, specfunc(chilled)),x);
     if not (bad = {}) then return FAIL; end if;

     # otherwise go ahead
     try
       c := solve({chill(b)},[x], 'useassumptions'=true) assuming op(as);
       postproc_for_solve((x,kb)->assert_deny_mb(x,pol,kb), warm(c), kb, as);
     catch "when calling '%1'. Received: 'side relations must be polynomials in (name or function) variables'":
       WARNING( sprintf( "siderels bug:\n\t'%s'\nwhen calling solve(%%1, %%2) assuming (%%3)"
                         , StringTools[FormatMessage](lastexception[2..-1])), b, x, as );
       return FAIL;
     end try;
   end proc;

   postproc_for_solve := proc(reduce_op, c, kb, as, $)
     local p, c0, c1;
     if c :: list and nops(c) = 0 then # false
       return FAIL;

     elif c :: list({relation, specfunc(`And`), `and`}) then # conjunction
       c0 := map(c -> if c::{specfunc(`And`),`and`} then op(c) else c end if,c);
       return foldr(((z,kb)->postproc_for_solve(reduce_op, z, kb, as)), kb, op(c0));

     elif c :: list then # disjunction
       if nops(c)=1 then
         return postproc_for_solve(reduce_op, op(1,c), kb, as);
       else
         return FAIL;
       end if;

     elif c :: {relation,specfunc(`Not`)} then # atom
       return reduce_op(c, kb);

     elif c :: specfunc(`piecewise`) then # try to make it into a conjunction
       p := Partition:-PWToPartition(c);
       p := apply(Partition:-Simpl:-remove_false_pieces,p) assuming op(as);
       c0, c1 := Partition:-Simpl:-single_nonzero_piece(p, _testzero=(x->x=[]));
       if not c0 :: identical(true) then
         if c1 :: relation then
         elif c1 :: list and nops(c1) = 1 then
             c1 := op(1,c1);
             if c1 :: list then c1 := op(c1) end if;
         else
             return FAIL;
         end if;
         try return postproc_for_solve(reduce_op, [ c0, c1 ], kb, as);
         catch "when calling '%1'. Received: 'cannot assume on a constant object'": NULL; end try;
       end if;
       return FAIL;
     end if;

     error "don't know what to do with %1 (in ctx %2, %3)", c, as, kb;
   end proc;

   # Simplify `bb' in context `as'
   simplify_in_context := proc(bb, as, $)
     local b;
     b := chill(bb);
     b := simplify(b) assuming op(as);
     warm(b);
   end proc;

   rel_coulditbe := proc(a,as0,$)
      local as := remove(x->x::`not`(`and`), as0);
      try
          coulditbe(a) assuming op(as);
      catch "when calling '%1'. Received: 'contradictory assumptions'" :
          # technically this means the KB was already contradictory, we
          # just didn't know?
          return false;
      catch "when calling '%1'. Received: 'side relations must be polynomials in (name or function) variables'":
          # This is seemingly a Maple bug - the condition could still be, but we
          # don't know, so conservatively return true.
          WARNING( sprintf( "siderels bug:\n\t'%s'\nwhen calling coulditbe(%%1) assuming (%%2)"
                          , StringTools[FormatMessage](lastexception[2..-1])), a, as0 );
          return true;
      catch "when calling '%3'. Received: 'when calling '%2'. Received: 'expression independent of, %0''":
          error expr_indp_errMsg(), a, as;
      catch "when calling '%2'. Received: 'expression independent of, %0'":
          error expr_indp_errMsg(), a, as;
      end try;
   end proc;

   expr_indp_errMsg := proc($)
       sprintf("Something strange happened(%s)\n"
               "\tin coulditbe(%%1) assuming %%2"
              ,StringTools[FormatMessage](lastexception[2..-1]))
   end proc;

   # Given a constraint "bb" on a KB "kb", this
   #   inserts either "bb" (if "pol" is true) or "Not bb" (otherwise)
   #   or, KB(Constrain(`if`(pol,bb,Not(bb))), kb)
   # Great deal of magic happens behind the scenes
   ModuleApply := proc(bb0::t_kb_atom, pol::identical(true,false), kb::t_kb)
    # Add `if`(pol,bb,Not(bb)) to kb and return the resulting KB.
    local as, bb, bbv, b, k, x, log_b, todo, kb0, ch;
    bb := bb0;

    if bb = pol then
      # Ignore literal true and Not(false).
      kb

    elif bb :: t_if_and_or_of(pol) then
      foldr(((b,kb) -> assert_deny_mb(b, pol, kb)), kb, op(bb))

    elif bb :: t_not then
      foldr(((b,kb) -> assert_deny_mb(b, not pol, kb)), kb, op(bb))

    else
      as := chill(kb_to_assumptions(kb, bb));
      bb := chill(bb);

      # try to evaluate under the assumptions, but some assumptions break
      # with eval, so remove any of those we tried to chill to prevent them breaking
      bb := subsindets(bb, relation, x->map(eval,x) assuming(op(as)));

      # Check that the new clause would not cause a contradictory
      # KB. If it does, then produce NotAKB.
      if not bad_assumption(bb) and not rel_coulditbe(`if`(pol,bb,not(bb)), as) then
          return NotAKB();
      end if;

      b := simplify_in_context(bb, as);

      # Look through kb for the innermost scope where b makes sense.
      k := select((k -> k :: Introduce(name, anything) and depends(b, op(1,k))),
                  kb);

      # If that scope is not precisely the trivial KB, then ..
      if nops(k) > 0 then
        x, k := op(op(1,k));
        # Found the innermost scope where b makes sense.
        # Reduce (in)equality between exp(A) and exp(B) to between A and B.
        b := try_improve_exp(b, x, as);

        # syntactic adjustment
        # If `b' is of a particular form (a bound on `x'), simplification
        # is in order
        if not b::`::` then
          todo := `if`(b :: t_bound_on(`x`), bound_simp, not_bound_simp);
          kb0 := todo(b,x,k,kb,pol,as);

          # If it succeeds, return that result
          if not kb0 :: identical(FAIL) then return kb0 end if;
        end if;
      end if;

      # Normalize `=` and `<>` constraints a bit.
      if not pol then
        b := negate_rel(b);
      end if;

      # If the name in the simple equality (if it is such) is not
      # on the lhs, then flip the equality
      if b :: t_constraint_flipped then b := (rhs(b)=lhs(b)) end if;

      # If `b' reduces to `true' in the KB environment then there is no need to
      # add it
      ch := chill(b);
      try
        if is(ch) assuming op(as) then return kb end if;
      catch "when calling '%1'. Received: 'side relations must be polynomials in (name or function) variables'":
        WARNING( sprintf( "siderels bug:\n\t'%s'\nwhen calling is(%%1) assuming (%%2)"
                        , StringTools[FormatMessage](lastexception[2..-1])), b, as );
      end try;

      # Add constraint to KB.
      KB(Constrain(b), op(kb))
    end if;
   end proc: # ModuleApply
  end module; # assert_deny

  # In order to hopefully produce a simplification,
  #   assert_deny will on some occasions repeatedly apply
  #   `simplify@ln` in the hopes of producing an 'improved'
  #   constraint. This metric gives the stopping condition
  # - when the simplification ceases to improve the constraint
  # - which is when the metric is made no less by the `simplify@ln`.
  # Since this is a `length`, this is strictly decreasing,
  #  so such a loop will 'provably' always terminate
  log_metric := proc(e, x, $)
    local m, L;
    m := select(depends, indets(e, 'exp(anything)'), x);
    length(subsindets(map2(op, 1, m), name, _->L));
  end proc;

  # This should be local to KB (or even to assert_deny) but it is used
  # by Domain.
  try_improve_exp := proc(b0, x, ctx, $)
        local b := b0, log_b;
        do
          try log_b := map(simplify@ln, b) assuming op(ctx); catch: break; end try;

          if log_metric(log_b, x) < log_metric(b, x)
             and (andmap(e->is(e,real)=true, log_b) assuming op(ctx)) then
            b := log_b;
          else
            break;
          end if;
        end do;
        b;
  end proc;

  # boolean_if should be equivalent to `if`, but it assumes
  # all its arguments are boolean conditions, so it basically
  # simplifies "cond and th or not cond and el"
  boolean_if := proc(cond, th, el, $)
    use
      a = ((x,y)-> `if`(x=true,y, `if`(x=false,x,
                   `if`(y=true,x, `if`(y=false,y, And(x,y)))))),
      o = ((x,y)-> `if`(x=false,y, `if`(x=true,x,
                   `if`(y=false,x, `if`(y=true,y, Or (x,y)))))),
      n = (x    -> `if`(x=false,true,
                   `if`(x=true,false,             Not(x))))
    in
      o(a(cond,th), a(n(cond),el))
    end use
  end proc;

  # Given that kb is an extension of kb0
  # (in that all the knowledge in kb0 is contained in kb)
  # then produces kb 'without' kb0.
  # Essentially this just applies coalesce_bounds
  # and then folds over the kb, at each step producing
  # a 'valid' KB by applying the smart constructor
  # (?) What is the role of of coalesce_bounds? why is it necessary?
  kb_subtract := proc(kb::t_kb, kb0::t_kb, $)
    local cut;
    cut := nops(kb) - nops(kb0);
    if cut < 0 or [op(cut+1..-1, kb)] <> [op(kb0)] then
      error "%1 is not an extension of %2", kb, kb0;
    end if;
    map(proc(k, $)
      local x, t;
      if k :: t_kb_Introduce then
        x, t := op(k);
        if t :: 'specfunc(AlmostEveryReal)' then
          t := [op(t), Bound(`>`, -infinity), Bound(`<`, infinity)];
          [genLebesgue, x,
           op([1,2], select(type, t, Bound(identical(`>`), anything))),
           op([1,2], select(type, t, Bound(identical(`<`), anything)))]
        else
          [genType, x, t]
        end if
      elif k :: t_kb_Let then
        [genLet, op(k)]
      elif k :: t_kb_Bound then
        [assert, op(2,k)(op(1,k),op(3,k))]
      elif k :: t_kb_Constrain then
        [assert, op(1,k)]
      end if
    end proc, [op(coalesce_bounds(KB(op(1..cut, kb))))])
  end proc;

  # checks that the given condition is fully contained in the given KB,
  # or that `kb => cond'
  kb_entails := proc(kb::t_kb, cond,$)
      local kb0;
      # nothing implies false (other than false, but kb is not false)
      if cond :: t_not_a_kb then
          false
      # if the condition is a kb, subtract from it directly
      elif cond :: t_kb then
          nops(kb_subtract(cond, kb)) = 0
      else # if the condition is an atom, check that asserting it has
           # no effect on the KB
          kb0 := assert(cond,kb);
          if kb0 :: t_kb then
              nops(kb_subtract(kb0,kb)) = 0
          else
              false
          end if;
      end if;
  end proc;

  # This essentially extracts all of the `Bound`s from a
  # KB and then re-inserts them 'directly' by applying their
  # knowledge to the rest of the KB. This may (?) produce
  # an invalid KB
  coalesce_bounds := proc(kb::t_kb, $)
    local lo, hi, eq, rest, k, x, t, b, s, r;

    # For every introduction in kb, remove bounds from
    # the introduction and store them seperately
    #  rest   maps variables to the stripped type
    #  lo,hi  map variables to low,high bounds
    for k in select(type, kb, t_intro) do
      x, t := op(k); # "x = t" is the intro

      # t := minus the lower,upper bounds in t
      b, t := selectremove(type, t, Bound(t_lo, anything));
      if nops(b) > 0 then lo[x] := op(1,b) end if;
      b, t := selectremove(type, t, Bound(t_hi, anything));
      if nops(b) > 0 then hi[x] := op(1,b) end if;
      rest[x] := [op(t)];
    end do;

    # Extract equality bounds, stored in eq (again a map from names)
    for k in select(type, kb, Bound(name, identical(`=`), anything)) do
      eq[op(1,k)] := op(3,k);
    end do;

    # Select `Bound`s in `Constrain`s, i.e. bounds on bound variables
    for k in select(type, kb, Bound(name, t_lo, anything)) do
      lo[op(1,k)] := subsop(1=NULL,k);
    end do;
    for k in select(type, kb, Bound(name, t_hi, anything)) do
      hi[op(1,k)] := subsop(1=NULL,k);
    end do;

    # 'coalesce' the bounds back into each clause of the KB
    map(proc(k, $)
      # when the clause is an introduction
      #  and there are equations on "x" which all evaluate to "a = a"
      #  remove that introduction and substitute
      #  bounds information about "x" back into things about "x"
      if k :: t_intro then
        x := op(1,k);
        if eq[x] = evaln(eq[x]) then
          subsop(2=op([2,0],k)(op(select(type, [lo[x],hi[x]], specfunc(Bound))),
                               op(rest[x])),
                 k);
        else
      #  otherwise, just replace the intro with a Let
          Let(x, eq[x]);
        end if

      # When the clause is a bound, erase it
      elif k :: t_kb_Bound
           and rest[op(1,k)] :: 'list' then
        NULL;

      # otherwise, do nothing
      else
        k;
      end if;
    end proc, kb);
  end proc;

  # The constraints do not include type assumptions
  kb_to_constraints := proc(kb::t_kb, $)::list(t_kb_atom);
    remove(type, kb_to_assumptions(kb), `::`);
  end proc;

  eval_kb := proc(e,kb::t_kb, $)
    foldl(eval, e, op(kb_to_equations(kb)));
  end proc;

  kb_assuming_mb := simpl -> proc(ee, kb::t_kb, failure, $)
    local e, as, e0;                                                         # for debugging

    e := eval_kb(ee,kb);                                                  `eval`;
    as := kb_to_assumptions(kb, e);
    e := chill(e);                                                        `chill`;
    as := chill(as);

    # The assumptions may be contradictory - I'm not sure why the right thing to
    # do is to return the original expression. Assuming false one can derive
    # anything - hence exception - so who calls this function under which
    # contexts that they expect `false` to mean something other than `false`?
    e0 := e;
    try e := simpl(e) assuming op(as); catch: e := failure(e0); end try;

    e := warm(e);                                            `warm (then expand@exp)`;
    eval(e, exp = expand @ exp);
  end proc;

  # Simplfies a given Hakaru term under knowledge of the
  # given KB. Does some magic to appease the Maple simplifier.
  # simplification might fail, in which case `failure(e)` where `e`
  # is the un-simplified (and chilled) expression is taken to be the result of
  # simplification. 'mb' for 'maybe'
  simplify_assuming_mb := kb_assuming_mb(simplify);

  simplify_assuming := proc(ee, kb::t_kb, $)
    simplify_assuming_mb(ee,kb,e->e);
  end proc;

  # like simplify_assuming but propgates the failure (as FAIL) instead of
  # silently consuming it and returning the unsimplified expression.  'f' for
  # 'failure'
  simplify_assuming_f := proc(ee, kb::t_kb, $)
    simplify_assuming_mb(ee,kb,e->FAIL);
  end proc;

  getType := proc(kb :: t_kb, x :: name, $)
    local res, over, bound;
    res := op([1,2], select(type, kb, 'Introduce(identical(x), anything)'));
    over := table([`<`=identical(`<`,`<=`), `<=`=identical(`<`,`<=`),
                   `>`=identical(`>`,`>=`), `>=`=identical(`>`,`>=`)]);
    for bound in select(type, kb, 'Bound(identical(x),
                                         identical(`<`,`<=`,`>`,`>=`),
                                         anything)') do
      res := remove(type, res, 'Bound'(over[op(2,bound)], 'anything'));
      res := op(0,res)(subsop(1=NULL,bound), op(res));
    end do;
    res
  end proc;

  # extract all introduced variables from a KB
  kb_to_variables := proc(kb :: t_kb, $)
    [op(map2(op, 1, select(type, kb, t_intro)))];
  end proc;

  bad_assumption_pw := (x->x::`=` and has(x,piecewise));
  # Returns true if the assumption is bad, false otherwise
  bad_assumption := proc(a, $)
    bad_assumption_pw(a) or
    # The case above is because the following takes forever:
    # simplify(piecewise(_a = docUpdate, aaa, bbb)) assuming i = piecewise(_a_
    # = docUpdate, zNew, idx[z, _a]), _a::integer, 0 <= _a, _a <= size[t]-1,
    # i::integer, 0 <= i, i <= size[as]-2, size[xs] = size[as]-1, size[z] =
    # size[t], docUpdate::integer, 0 <= docUpdate, docUpdate <= size[z]-1
    ( a :: `=` and
      ormap(f->f(a)::name,[lhs,rhs]) and
      indets(a,'{specindex,specfunc}'(chilled))<>{} )
    # These are dealt with otherwise and aren't understood by Maple
  end proc;

  # Note that this returns potentially any number of operands.
  kb_atom_to_assumptions := proc(k, $)
    local x, x_ty, x_typ;
    if k :: t_intro then
      x, x_ty := op(k);
      x_typ := htype_to_property(x_ty);
      (x :: x_typ),
       op(`if`(x_typ<>TopProp, map((b -> `if`(nops(b)>=2,[op(1,b)(x, op(2,b))],[])[]), x_ty),[]))
    elif k :: t_kb_Let then
      `=`(op(k))
    elif k :: t_kb_Bound then
      op(2,k)(op(1,k), op(3,k))
    elif k :: t_kb_Constrain then
      op(1,k)
    else
      NULL # Maple doesn't understand our other types
    end if
  end proc;

  kb_to_assumptions := proc(kb, e:={}, to_remove := bad_assumption_pw, $)
    local n, as;
    as := remove(to_remove,
          map( kb_atom_to_assumptions ,
             [op(coalesce_bounds(kb))
             ,seq(Constrain(n::nonnegint), n in indets(e, 'specfunc(size)'))]));
  end proc;

  # extract Lets and equality constraints (only!) from a KB
  kb_to_equations := proc(kb, $)
    local lets, constraints;

    lets := map2(subsop, 0=`=`, [op(select(type, kb, t_kb_Let))]);
    constraints := map(op, select(type, kb, 'Constrain(anything = anything)'));
    [op(lets), op(constraints)]
  end proc;

  htype_to_property := proc(t::t_type, $)
    if t :: 'specfunc({AlmostEveryReal, HReal})' then real
    elif t :: 'specfunc({HInt})' then integer
    else TopProp end if
  end proc;

  # See kb_Partition
  kb_piecewise := proc(e :: specfunc(piecewise), kb :: t_kb, doIf, doThen)
    Partition:-PartitionToPW(
        kb_Partition(Partition:-PWToPartition(e, _rest), kb, doIf, doThen)
        ) ;
  end proc;

  # A sort of map over a Partition with the given KB as a context, such that:
  #    kb_Partition( PARTITION ( Piece( c_i , v_i ), .. )  , kb, doIf, doThen )
  #    =
  #    PARTITION ( Piece( doIf(c_i, kb), doThen(v_i, assert(c_i, kb)) ) )
  # Semantics originally given here:
  #  https://github.com/hakaru-dev/hakaru/commit/6f1c1ea2d039a91c157462f09f15760c98884303
  kb_Partition:= proc(e::Partition, kb::t_kb, doIf, doThen,$)::Partition;
  local br;
    #Unlike `piecewise`, the conditions in a Partition are necessarily
    #disjoint, so the `update` used in kb_piecewise isn't needed. We may
    #simply `assert` the condition (i.e., roll it into the kb) without
    #needing to `assert` the negation of all previous conditions.

    Partition:-Amap([(x,_) -> doIf(x, kb), (x,c) -> doThen(x, assert(c,kb)), z -> z], e);
  end proc;

  # Like convert(e, 'list', `*`) but tries to keep the elements positive
  list_of_mul := proc(e, $)
    local rest, should_negate, can_negate, fsn;
    rest := convert(e, 'list', `*`);
    rest := map((f -> [f, signum(f),
                       `if`(f::'{specfunc({Sum,sum}),anything^odd}',
                            applyop(`-`,1,f),
                       `if`(f::'specfunc(ln)',
                            applyop(`/`,1,f),
                            -f))]),
                rest);
    should_negate, rest := selectremove(type, rest,
      '[anything, -1, Not(`*`)]');
    if nops(should_negate) :: even then
      [seq(op(3,fsn), fsn=should_negate),
       seq(op(1,fsn), fsn=rest)]
    else
      can_negate, rest := selectremove(type, rest,
        '[{`+`, specfunc({Sum,sum,ln}), `^`}, Not(1), Not(`*`)]');
      if nops(can_negate) > 0 then
        [seq(op(3,fsn), fsn=should_negate),
         op([1,3], can_negate),
         seq(op(1,fsn), fsn=subsop(1=NULL, can_negate)),
         seq(op(1,fsn), fsn=rest)]
      else
        [seq(op(3,fsn), fsn=subsop(-1=NULL, should_negate)),
         op([-1,1], should_negate),
         seq(op(1,fsn), fsn=rest)]
      end if
    end if
  end proc;

  # Given an expression containing products and sums, i.e. polynomials
  # and a function , applies this expression to each factor and summand
  for_poly := proc(e, f, $)
    if e :: '{`+`,`*`}' then map(for_poly, e, f)
    elif e :: 'specfunc({product,Product,sum,Sum})' then
      applyop(for_poly, 1, e, f)
    else f(e)
    end if
  end proc;

  # Computes the range of (possible values of) a Hakaru Int,
  # given a Hakaru type for that Int
  range_of_HInt := proc(t :: And(specfunc(HInt), t_type), $)
       op(1, [op(map((b -> `if`(op(1,b)=`>`, floor(op(2,b))+1, op(2,b))),
                     select(type, t, Bound(t_lo,anything)))),
              -infinity])
    .. op(1, [op(map((b -> `if`(op(1,b)=`<`, ceil(op(2,b))-1, op(2,b))),
                     select(type, t, Bound(t_hi,anything)))),
              infinity]);
  end proc;

  # Avoid FAILure modes of the assumption system
  # This transforms {size,idx}(a,b,..x) to {size,id}[a,b,..x] (that is, a
  # subscript operator, which doesn't evaluate) and back. Some functions
  # evaluating causes simplification to fail because information is lost
  chilled := '{size, idx}';

  chill := e -> subsindets(e, specfunc(chilled), c->op(0,c)[op(c)]);
  warm := e -> subsindets(e, specindex(chilled), c->map(warm, op(0,c)(op(c))));

  # The KB constructors are local, but sometimes for debugging purposes one
  # would like to construct the KB directly. This converts the global names
  # Introduce,Constrain,Let, and Bound to KB:-<..> and replaces the 0-th operand
  # with KB:-KB.
  build_unsafely := proc(kb,$)
    KB(op(subs([ :-Introduce=Introduce, :-Let=Let, :-Bound=Bound, :-Constrain=Constrain ], kb)));
  end proc;

  ModuleLoad := proc($)
    ModuleUnload();
    Hakaru; # Make sure the KB module is loaded, for the type t_type
    TypeTools[AddType](t_kb,
      'specfunc({
         Introduce(name, t_type),
         Let(name, anything),
         Bound(name, identical(`<`,`<=`,`>`,`>=`,`=`), anything),
         Constrain({`::`, boolean, `in`, specfunc(anything,{Or,Not})})
       }, KB)');

    # KB 'atoms' , i.e. single pieces of knowledge, in "maple form".
    # Note that boolean already includes `Bound`s in the form `x R y`
    TypeTools[AddType](t_kb_atom,
      '{`::`, boolean, `in`, specfunc(anything,{Or,Not,And})}');
    TypeTools[AddType](t_kb_atoms,list(t_kb_atom));

    # The 'false' KB, produced when a contradiction arises in a KB
    TypeTools[AddType](t_not_a_kb, 'specfunc(NotAKB)');

    # Something that might be a KB, or is the false KB
    TypeTools[AddType](t_kb_mb, '{t_kb, t_not_a_kb}');

    # Prevent expand(product(f(i),i=0..n-1))
    # from producing (product(f(i),i=0..n))/f(n)
    `expand/product` := overload([
      proc(ff, rr::name=And(`+`,Not(`+`(Not(integer))))..anything)
        option overload(callseq_only);
        thaw(`expand/product`(ff, applyop(freeze, [2,1], rr)))
      end proc,
      proc(ff, rr::name=anything..And(`+`,Not(`+`(Not(integer)))))
        option overload(callseq_only);
        thaw(`expand/product`(ff, applyop(freeze, [2,2], rr)))
      end proc,
      `expand/product`]);

    # Prevent simplify(product(1-x[i],i=0..n-1))
    # from producing (-1)^n*product(-1+x[i],i=0..n-1)
    `simplify/int/simplify` := overload([
      proc (f::identical(product), a::`+`, r::{name,name=anything}, $)
        option overload(callseq_only);
        return 'f'(a,r);
      end proc,
      `simplify/int/simplify`]);

    # Prevent product(1/(n-i),i=0..n-1) and product(1/(i-n),n=0..i-1)
    # from evaluating to 0 (?!)
    # by preventing `product/indef/indef`(n-i,i)
    # from evaluating to (-1)^i * GAMMA(i-n)
    `product/indef/indef` := overload([
      proc(expr, i, $)
        option overload;
        local a, b;
        if expr :: `*` then
          `*`(op(map(`product/indef/indef`, list_of_mul(expr), i)))
        elif expr :: 'anything^freeof(i)' then
          `product/indef/indef`(op(1,expr),i)^op(2,expr)
        elif has(expr,i) and sign(expr,i) :: negative then
          if ispoly(expr, linear, i, 'b', 'a') then
            (-a)^i / GAMMA(1-i-b/a)
          else
            error FAIL
          end if
        else
          error "invalid input: cannot product/indef/indef(%1, %2)", expr, i
        end if
      end proc,
      `product/indef/indef`]);

    # Convert GAMMA(x)*GAMMA(y)/GAMMA(x+y) to Beta(x,y)
    `convert/Beta` := proc(e, $)
      subsindets(e, 'And(`*`, Not(`*`(Not(specfunc(GAMMA)))),
                              Not(`*`(Not(1/specfunc(GAMMA)))))',
        proc(p, $)
          local s, t, r, i, j, x, y;
          s, r := selectremove(type, convert(p,'list',`*`), 'specfunc(GAMMA)');
          t := map2(op, [1,1], select(type, {op(r)}, '1/specfunc(GAMMA)'));
          for i from nops(s) to 2 by -1 do
            x := op([i,1],s);
            for j from i-1 to 1 by -1 do
              y := op([j,1],s);
              if x + y in t then
                s := subsop(j=GAMMA(x+y), i=Beta(x,y), s);
                break;
              end if;
            end do;
          end do;
          `*`(op(s), op(r));
        end proc);
    end proc;
  end proc;

  ModuleUnload := proc($)
    map(proc(x::uneval) try eval(x) catch: NULL; end try end proc,
         ['TypeTools[RemoveType](t_kb)']);
  end proc;
  ModuleLoad();
end module; # KB

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
process ACV_setupScoreboards {
/scoreboard objectives add ACV_COLOR dummy
/scoreboard objectives add ACV_DIRECTION dummy
/scoreboard objectives add ACV_CONSTANT dummy
/scoreboard players set 2 ACV_CONSTANT 2
/scoreboard players set 5 ACV_CONSTANT 5
/scoreboard players set 4 ACV_CONSTANT 4
/scoreboard objectives add ACV_TP_FROM dummy
/scoreboard objectives add ACV_TP_UNEVEN dummy
/scoreboard objectives add ACV_TP_TO dummy
// Rotation of the Portal the entity goes in
/scoreboard objectives add ACV_TP_IR dummy
// Rotation of the Portal the entity comes out
/scoreboard objectives add ACV_TP_OR dummy
// Rotation Difference of IR and OR
/scoreboard objectives add ACV_TP_RD dummy
/scoreboard objectives add ACV_CUBE_ID dummy
/scoreboard objectives add OnGround dummy
}


text:
# File: MapleCodePrime.mpl
#
# Description:
#  This file contains code written by external developers -- users of mapleprimes.com
# Author:
#  Alec, Mihailovs
#
#
# Date:
#  27/12/2016
#
# License:
#  Contact the authors.
#
# Copyright (c) Alec Mihailovs
# All rights reserved.

RigidMotionsMaplePrimesCode := module()
  option package;

  export Isort;

# Procedure: Isort
#   Sort elements and return their indices in original order after sorting
#
# Author:
#   Alec Mihailovs - http://www.mapleprimes.com/posts/43507-Sorting-With-Indices#comment80473
# Parameters:
#   L            - a container to be sorted.
#
# Output:
#   indices of original order after sorting, sorted collection.
Isort:=proc(L)
  local a;
  a := sort([$1..nops(L)],(i,j)->L[i]<=L[j]);
  return a, [seq(L[i],i=a)];
end:

end module;

max_stars_repo_name: copyme/MapleTools

max_issues_repo_name: copyme/MapleTools

max_forks_repo_name: copyme/MapleTools


text:
print "Hello MiniPL\n";
assert(false);

max_stars_repo_name: nicohi/mini-pl-interpreter

max_issues_repo_name: nicohi/mini-pl-interpreter

max_forks_repo_name: nicohi/mini-pl-interpreter


text:
with(RegularChains):
with(SemiAlgebraicSetTools):
with(ListTools):

EvalCertification := module()
option package;

export HolderInformationForPolynomial,
	   HolderInformationForExponential,
	   HolderInformationForRationalPolynomial,
	   EstimateRootsAndCertifyEvaluations,
	   ExpandToImaginaryAndReal;

#ExpandToImaginaryAndReal,
local SqrFreeFactorization,
	  SingleEqualityToEstimate,
	  BoxEqualitiesToEstimate,
	  EstimateAndClassifySolutions,
	  HolderDataToCertificationTolerance,
	  EvaluateFunctionAtRoot;

# Compute a square free factorization of a polynomial
# along with multiplicities
SqrFreeFactorization := proc(F)
	local FsqrFree, fa, CurrOutput, j, k; FsqrFree := sqrfree(F);
	FsqrFree := FsqrFree[2]; CurrOutput := [];
	for j to nops(FsqrFree) do
		fa := ifelse(type(factor(FsqrFree[j][1]), `+`), [FsqrFree[j][1]], convert(factor(FsqrFree[j][1]), list));
		for k to nops(fa) do
			CurrOutput := [op(CurrOutput), [fa[k], FsqrFree[j][2]]];
		end do;
	end do;
	return CurrOutput;
end proc;

# Expand a list of polynomials to list with real and
# imaginary parts separated
ExpandToImaginaryAndReal := proc(F)
	local vars, imaginary_variables, expanded_equations, i;
	vars := indets(F);
	imaginary_variables := [seq(b[i], i = 1 .. nops(vars))];
	expanded_equations := expand(subs([seq(vars[i] = Complex(vars[i], imaginary_variables[i]), i = 1 .. nops(vars))], F));
	return [[op(subs(I = 0, expanded_equations)), op([seq(coeff(expanded_equations[i], I), i = 1 .. nops(expanded_equations))])], [seq(vars[i]^2 + imaginary_variables[i]^2 - 1, i = 1 .. nops(vars))]];
end proc;

SingleEqualityToEstimate := proc(equality)
	return ifelse(type(rhs(equality), 'list'), 1/2*rhs(equality)[2] + 1/2*rhs(equality)[1], rhs(equality));
end proc;

BoxEqualitiesToEstimate := proc(equalities)
	local reduced_equalities;
	reduced_equalities := map(SingleEqualityToEstimate, equalities);
	return Complex(reduced_equalities[1], reduced_equalities[2]);
end proc;

# Estimate the solution of univarite polynomial F and classify
# by whether the imaginary norm is > 1, < 1, or = 1
EstimateAndClassifySolutions := proc(F, epsilon)
	# INPUT: Squarefree polynomial F in format [polynomial,power],
	# estimation error epsilon

	# OUTPUT: [list_of_estimates_for_norm_1_solutions,list_of_estimates_for_norm_<_1_solutions,
	# 		   list_of_estimates_for_norm_>_1_solutions]
	local actual_polynomial, system_with_imaginary_parts, equations,
		  norms, R, norm_1_boxes, norm_1_values, norm_greater_boxes,
		  norm_greater_values, norm_less_boxes, norm_less_values, i,
		  norm_1_values_refined, norm_greater_values_refined, norm_less_values_refined;
	actual_polynomial := [F[1]];
	system_with_imaginary_parts := ExpandToImaginaryAndReal(actual_polynomial);
	equations := system_with_imaginary_parts[1];
	norms := system_with_imaginary_parts[2];
	R := PolynomialRing([seq(indets(norms)[i], i = 1 .. nops(indets(norms)))]);
	norm_1_boxes := RealRootIsolate([op(equations)], [], [], [], R, 'abserr' = epsilon);
	norm_1_values := map(BoxValues, norm_1_boxes, R);
	norm_1_values_refined := map(BoxEqualitiesToEstimate, norm_1_values);
	norm_1_values_refined := [seq(Record("value"=norm_1_values_refined[i]), i = 1 .. nops(norm_1_values_refined))];
	return [op(norm_1_values_refined)];
end proc;

# Compute local Lipschitz constant for univarite polynomial F at a point
# given an estimate for that point and a bound on that estimate's accuracy
HolderInformationForPolynomial := proc(F, PointEstimate, EstimateAccuracy)
	# INPUT: Univariate polynomial F,
	# an estimate PointEstimate of a point in F's domain (could be real or complex),
	# an upper bound EstimateAccuracy on the distance from PointEstimate to the
	# actual root

	# OUTPUT: A Record of the form [exponent,constant] such that
	# F is certifiably Holder with those values
	local deg, Derivatives, HolderConstant, HolderExponent, i;
	deg := degree(F);
	if deg > 0 then
		Derivatives := Vector();
		Derivatives(1) := diff(F, indets(F)[1] $ 1);
		for i from 2 to deg do
	    	Derivatives(i) := diff(Derivatives(i - 1), indets(F)[1] $ 1);
		end do;
		Derivatives := convert(Derivatives, list);
		HolderConstant := add((EstimateAccuracy)^(i - 1)*abs(subs(indets(F)[1] = PointEstimate, Derivatives[i]))/(i - 1)!, i = 1 .. deg);
	else HolderConstant := 0;
	end if;
	HolderExponent := 1;
	return Record('exponent'=HolderExponent,'constant'=HolderConstant,'avoid_roots'=false,'max_degree'=deg-1);
end proc;

# Holder information for exponentials of the form x^alpha
HolderInformationForExponential := proc(function_exponent)
	local output_function;
	output_function := proc(F, PointEstimate, EstimateAccuracy)
		return Record('exponent'=function_exponent,'constant'=1,'avoid_roots'=false,'max_degree'=1);
	end proc;
	return output_function;
end proc;

HolderInformationForRationalPolynomial := proc(P, PointEstimate, EstimateAccuracy,domain_estimate:=false)
    # INPUT: Univariate rational polynomial  of the form F/G

    # OUTPUT: A record of the form [exponent, constant] such that
    # F/G is certifiably Holder with those values
	local F, G, differential_numerator, numerator_constant, max_value,
		  denominator_constant, min_value,FactorsAndMultiplicities,ClassifiedRoots,MINIMUM_NONZERO;
	if domain_estimate then
		FactorsAndMultiplicities := SqrFreeFactorization(denom(simplify(P)));
		ClassifiedRoots := map(EstimateAndClassifySolutions, FactorsAndMultiplicities, EstimateAccuracy);
		ClassifiedRoots := FlattenOnce(ClassifiedRoots);
		return Record('avoid_roots'=[seq(ClassifiedRoots[i]:-value,i=1..nops(ClassifiedRoots))]);
	end if;

	MINIMUM_NONZERO := 3*10^(-10);
	F := numer(simplify(P));
	G := denom(simplify(P));

	if type(F, integer) then
		differential_numerator := -diff(G, indets(G)[1] $ 1)*F;
	elif type(G, integer) then
		differential_numerator := G*diff(F, indets(F)[1] $ 1);
	else differential_numerator := -diff(G, indets(G)[1] $ 1)*F + G*diff(F, indets(F)[1] $ 1);
	end if;
	numerator_constant := HolderInformationForPolynomial(differential_numerator, PointEstimate, EstimateAccuracy):-constant;

	if type(differential_numerator, integer) then
		max_value := abs(differential_numerator) + numerator_constant*EstimateAccuracy;
	else max_value := abs(subs(indets(differential_numerator)[1] = PointEstimate, differential_numerator)) + numerator_constant*EstimateAccuracy;
	end if;
	if type(G, integer) then
		denominator_constant := 1;
	else denominator_constant := HolderInformationForPolynomial(G^2, PointEstimate, EstimateAccuracy):-constant;
	end if;
	if evalb(numerator_constant = denominator_constant) then
		return MINIMUM_NONZERO;
	end if;
	if type(G, integer) then
		min_value := G^2 - denominator_constant*EstimateAccuracy;
	else min_value := abs(subs(indets(G)[1] = PointEstimate, G))^2 - denominator_constant*EstimateAccuracy;
	end if;

	return Record('exponent'=1,'constant'=max_value/min_value,'max_degree'=1);
end proc;

HolderDataToCertificationTolerance := proc(HolderData, ErrorTolerance)
	local HolderConstant, HolderExponent,RoundedDigitsTolerance,RoundedCertNegativeExponent;
	HolderExponent := HolderData:-exponent;
	HolderConstant := HolderData:-constant;
	RoundedDigitsTolerance := max(-floor(log[10](ErrorTolerance^HolderData:-max_degree)),1) + 1;
	if evalf[RoundedDigitsTolerance](HolderConstant) <= 2*10^(-RoundedDigitsTolerance) then
		return 0;
	end if;
	RoundedCertNegativeExponent := floor(-log[10](ErrorTolerance/HolderConstant)/HolderExponent);
	return min(1/10^RoundedCertNegativeExponent, ErrorTolerance, 1);
end proc;

EvaluateFunctionAtRoot := proc(Function, Root, Precision)
	local SymbolicEval,DigitsBeforeDecimalPoint,DigitsAfterDecimalPoint;
	if type(Function,procedure) then
		SymbolicEval := Function(Root);
	else
		SymbolicEval := subs(indets(Function)[1]=Root,Function);
	end if;
	if evalb(SymbolicEval=0) then
		return 0;
	end if;
	# Precision calculation: Figure out necessary number of digits
	# in front of the decimal point, figure out the same for after
	# the decimal point, then add them together.
	DigitsBeforeDecimalPoint := max(ceil(evalf[1](log[10](abs(SymbolicEval)))),1);
	DigitsAfterDecimalPoint := max(-floor(log[10](Precision)),0);
	return evalf[max(1,DigitsBeforeDecimalPoint+DigitsAfterDecimalPoint+1)](SymbolicEval);
end proc;

EstimateRootsAndCertifyEvaluations := proc(PolynomialToSolve, FunctionsToEvaluate, HolderEstimationProcedure, ErrorTolerance, AbandonThreshold := 10^(-100))
	local FactorsAndMultiplicities, Multiplicities, ClassifiedRoots,
	      RootsWithoutAdditionalInformation, HolderBounds, HolderCertificationBounds,
	      SharpenedRootTolerance, ClassifiedRootsWithMultiplicities, WorkingTolerance, ActualTolerance,
	      i, j,arguments_for_output,root_information_for_output,evaluation_information_for_output,current_evals,
	      minimum_gap,current_roots,holder_function,function_roots,estimated_gap_values;
	FactorsAndMultiplicities := SqrFreeFactorization(expand(PolynomialToSolve));
	Multiplicities := [seq(FactorsAndMultiplicities[i][2], i = 1 .. nops(FactorsAndMultiplicities))];
	SharpenedRootTolerance := 0;
	# We need to make sure all estimates are at least 1/10th sharper
	# than the desired number because rounding error in the last reported digit
	# could result in that amount of error
	ActualTolerance := ErrorTolerance - (1/10)*ErrorTolerance;
	WorkingTolerance := max(1/10,ActualTolerance);

	# First Holder information step: Tightening the precision
	# to avoid places where the evaluation functions aren't
	# defined. We need that there are no undefined points
	# with 2*WorkingTolerance of the function to solve's
	# estimated roots
	for i to nops(FunctionsToEvaluate) do
		minimum_gap := 0;
		if type(HolderEstimationProcedure,list) then
			holder_function := HolderEstimationProcedure[i];
		else
			holder_function := HolderEstimationProcedure;
		end if;
		# Note: Minimum gap often contains square roots and Maple doesn't
		# automatically compare expressions with abs in them symbolically.
		# For minimum_gap we can just square, for others we need to evaluate
		# and be careful with precision
		while minimum_gap^2 <= (6*WorkingTolerance)^2 and WorkingTolerance > AbandonThreshold do
			if minimum_gap^2 > 0 then
				WorkingTolerance := min(WorkingTolerance/2,minimum_gap/(6*WorkingTolerance));
			end if;
			function_roots := holder_function(FunctionsToEvaluate[i],1.1,WorkingTolerance,true):-avoid_roots;
			if evalb(function_roots=false) then
				break;
			end if;
			ClassifiedRoots := map(EstimateAndClassifySolutions, FactorsAndMultiplicities, WorkingTolerance);
			ClassifiedRoots := FlattenOnce(ClassifiedRoots);
			current_roots := [seq(ClassifiedRoots[j]:-value, j = 1 .. nops(ClassifiedRoots))];
			estimated_gap_values := FlattenOnce([ seq([seq(abs(function_roots[j] - current_roots[k]),k=1..nops(current_roots))],j=1..nops(function_roots)) ]);
			minimum_gap := min(estimated_gap_values);
			# Need to round down to some rational representation
			# to keep things rational
			minimum_gap := convert(evalf[10](minimum_gap)-2*10^(-10),rational);
		end do
	end do;

	# At this point we've reduced to a good enough tolerance to
	# find local Holder information while dodging points not in the domains
	# of the functions to evaluate. Proceed to find that information, sharpen tolerance,
	# and give final results.
	while SharpenedRootTolerance = 0 and AbandonThreshold < WorkingTolerance do
		ClassifiedRoots := map(EstimateAndClassifySolutions, FactorsAndMultiplicities, WorkingTolerance);
		ClassifiedRoots := FlattenOnce(ClassifiedRoots);
		RootsWithoutAdditionalInformation := [seq(ClassifiedRoots[i]:-value, i = 1 .. nops(ClassifiedRoots))];
		if type(HolderEstimationProcedure, list) then
			if nops(HolderEstimationProcedure) <> nops(FunctionsToEvaluate) then
				error "Number of Holder estimation procedures is different than the number of functions to evaluate.";
			end if;
			HolderBounds := [seq(map[2](HolderEstimationProcedure[i], FunctionsToEvaluate[i], RootsWithoutAdditionalInformation, 2*WorkingTolerance), i = 1 .. nops(FunctionsToEvaluate))];
		else
			HolderBounds := [seq(map[2](HolderEstimationProcedure, FunctionsToEvaluate[i], RootsWithoutAdditionalInformation, 2*WorkingTolerance), i = 1 .. nops(FunctionsToEvaluate))];
		end if;
		HolderBounds := FlattenOnce(HolderBounds);
		HolderCertificationBounds := map(HolderDataToCertificationTolerance, HolderBounds, min(ActualTolerance,WorkingTolerance/2));
		SharpenedRootTolerance := min(HolderCertificationBounds);
		WorkingTolerance := 1/10*WorkingTolerance;
	end do;
	if SharpenedRootTolerance = 0 then
		error "Insufficient precision available to certify evaluations.";
	end if;
	ClassifiedRoots := map(EstimateAndClassifySolutions, FactorsAndMultiplicities, SharpenedRootTolerance);
	ClassifiedRoots := FlattenOnce(ClassifiedRoots);
	arguments_for_output := Record("polynomial"=PolynomialToSolve,"evaluated_functions"=FunctionsToEvaluate,"error_tolerance"=ActualTolerance);
	root_information_for_output := Record("root_values"=[seq(ClassifiedRoots[i]:-value, i=1..nops(ClassifiedRoots))], "root_multiplicities"=Multiplicities);
	evaluation_information_for_output := Record();
	for i to nops(FunctionsToEvaluate) do
		current_evals := Record("evaluations_functions_"||i=map[2](EvaluateFunctionAtRoot,FunctionsToEvaluate[i],root_information_for_output:-root_values,ActualTolerance));
		evaluation_information_for_output := Record[evaluation_information_for_output,current_evals]();
	end do;
	return Record[arguments_for_output,root_information_for_output,evaluation_information_for_output]();
end proc;

end module:

max_stars_repo_name: P-Edwards/EvalCertification

max_issues_repo_name: P-Edwards/EvalCertification

max_forks_repo_name: P-Edwards/EvalCertification


text:
var $const1 <[13] i8> = [ 104, 101, 108, 108, 111, 32, 119, 111, 114, 108, 100, 10, 0 ]
func $printf (var %p1 <* i8>)void
func $main ( ) i32 {
   call &printf (addrof a32 $const1)
   return (constval i32 0) }

 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)
(* prefix:
  gga_k_ol2_params *params;

  assert(p->params != NULL);
  params = (gga_k_ol2_params * )(p->params);
*)

ol2_f := x ->
  + params_a_aa
  + params_a_bb*x^2/72.0
  + params_a_cc*x/(2^(1/3) + 4*x):

f := (rs, zeta, xt, xs0, xs1) -> gga_kinetic(ol2_f, rs, zeta, xs0, xs1):


text:
program arraytest2;
var a : array[1] of integer;
var b : array[0] of char;

begin
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing

max_forks_repo_name: Hiroya-W/lang-processing


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)

$define gga_c_pbe_params
$include "gga_c_pbe.mpl"

(* in the paper we have beta_a = 0.066725 *)
beta_a := 0.066724550603149220:
beta_b := 0.1:
beta_c := 0.1778:

(* we redefine beta here *)
(* this is the Hu and Langreth expression *)
mbeta := (rs, t) -> beta_a*(1 + beta_b*rs)/(1 + beta_c*rs):


text:
program assigninteger	{sample21}
var n: integer;
begin
 n := 1
end.

max_issues_repo_name: shouth/LanguageProcessing


text:
program sample28p;
procedure p;
    begin
        writeln('Hello!')
    end;
procedure q;
    begin
        writeln('Everyone!')
    end;
begin
    call p;
    call q
end.


text:
repeat process ACV_smoothTeleportation {
/scoreboard players reset @e ACV_TP_FROM
/scoreboard players reset @e ACV_TP_UNEVEN
/scoreboard players reset @e ACV_TP_TO
/scoreboard players reset @e ACV_TP_IR
/scoreboard players reset @e ACV_TP_OR
/scoreboard players reset @e ACV_TP_RD

/execute @e[type=area_effect_cloud,name=ACV_PortalIn] ~ ~ ~ scoreboard players operation @e[tag=!ACV_TpProof,dy=0] ACV_TP_FROM = @e[type=area_effect_cloud,name=ACV_PortalIn,r=0,c=1] ACV_COLOR

// Calulate the target portal
/execute @e[score_ACV_TP_FROM_min=0] ~ ~ ~ scoreboard players operation @e[score_ACV_TP_FROM_min=0,r=0,c=1] ACV_TP_TO = @e[score_ACV_TP_FROM_min=0,r=0,c=1] ACV_TP_FROM
/execute @e[score_ACV_TP_FROM_min=0] ~ ~ ~ scoreboard players operation @e[score_ACV_TP_FROM_min=0,r=0,c=1] ACV_TP_UNEVEN = @e[score_ACV_TP_FROM_min=0,r=0,c=1] ACV_TP_FROM
/scoreboard players operation @e[score_ACV_TP_FROM_min=0] ACV_TP_UNEVEN %= 2 ACV_CONSTANT
/scoreboard players add @e[score_ACV_TP_UNEVEN_min=0,score_ACV_TP_UNEVEN=0] ACV_TP_TO 1
/scoreboard players remove @e[score_ACV_TP_UNEVEN_min=1,score_ACV_TP_UNEVEN=1] ACV_TP_TO 1

// Assign In-Rotation
/execute @e[type=area_effect_cloud,name=ACV_PortalIn] ~ ~ ~ scoreboard players operation @e[score_ACV_TP_TO_min=0,dy=0] ACV_TP_IR = @e[type=area_effect_cloud,name=ACV_PortalIn,r=0,c=1] ACV_DIRECTION

// Lightbridges must be rotated to look away from the portal, because the idle bridge is a fresh summon that always looks south
/scoreboard players set @e[tag=ACV_RotationProof] ACV_TP_IR -1

// Assign Out-Rotation

// Bug Workaround: https://bugs.mojang.com/browse/MC-89667
// /tp @e[score_ACV_TP_TO_min=0,score_ACV_TP_TO=0] @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=0,score_ACV_COLOR=0]
// /tp @e[score_ACV_TP_TO_min=1,score_ACV_TP_TO=1] @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=1,score_ACV_COLOR=1]
// /tp @e[score_ACV_TP_TO_min=2,score_ACV_TP_TO=2] @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=2,score_ACV_COLOR=2]
// /tp @e[score_ACV_TP_TO_min=3,score_ACV_TP_TO=3] @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=3,score_ACV_COLOR=3]
///execute @e[name=ACV_PortalOut] ~ ~ ~ scoreboard players operation @e[score_ACV_TP_TO_min=0,tag=!ACV_TpProof,dy=0] ACV_TP_OR = @e[r=0,c=1] ACV_DIRECTION

// Bug Workaround: https://bugs.mojang.com/browse/MC-89667
/scoreboard players operation @e[score_ACV_TP_TO_min=0,score_ACV_TP_TO=0] ACV_TP_OR = @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=0,score_ACV_COLOR=0] ACV_DIRECTION
/scoreboard players operation @e[score_ACV_TP_TO_min=1,score_ACV_TP_TO=1] ACV_TP_OR = @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=1,score_ACV_COLOR=1] ACV_DIRECTION
/scoreboard players operation @e[score_ACV_TP_TO_min=2,score_ACV_TP_TO=2] ACV_TP_OR = @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=2,score_ACV_COLOR=2] ACV_DIRECTION
/scoreboard players operation @e[score_ACV_TP_TO_min=3,score_ACV_TP_TO=3] ACV_TP_OR = @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=3,score_ACV_COLOR=3] ACV_DIRECTION

// Rotation Calculation
/execute @e[score_ACV_TP_IR_min=0,score_ACV_TP_OR_min=0] ~ ~ ~ scoreboard players operation @e[score_ACV_TP_IR_min=0,score_ACV_TP_OR_min=0,r=0,c=1] ACV_TP_RD = @e[score_ACV_TP_IR_min=0,score_ACV_TP_OR_min=0,r=0,c=1] ACV_TP_IR
/execute @e[score_ACV_TP_IR_min=0,score_ACV_TP_OR_min=0] ~ ~ ~ scoreboard players operation @e[score_ACV_TP_IR_min=0,score_ACV_TP_OR_min=0,r=0,c=1] ACV_TP_RD -= @e[score_ACV_TP_IR_min=0,score_ACV_TP_OR_min=0,r=0,c=1] ACV_TP_OR
/execute @e[score_ACV_TP_IR_min=0,score_ACV_TP_OR_min=0] ~ ~ ~ scoreboard players add @e[score_ACV_TP_IR_min=0,score_ACV_TP_OR_min=0,r=0,c=1] ACV_TP_RD 4
/execute @e[score_ACV_TP_IR_min=0,score_ACV_TP_OR_min=0] ~ ~ ~ scoreboard players operation @e[score_ACV_TP_IR_min=0,score_ACV_TP_OR_min=0,r=0,c=1] ACV_TP_RD %= 4 ACV_CONSTANT

// Teleport the entity to the other portal

// Bug Workaround: https://bugs.mojang.com/browse/MC-89667
/execute @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=0,score_ACV_COLOR=0] ~ ~ ~ teleport @e[score_ACV_TP_TO_min=0,score_ACV_TP_TO=0] ~ ~ ~
/execute @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=1,score_ACV_COLOR=1] ~ ~ ~ teleport @e[score_ACV_TP_TO_min=1,score_ACV_TP_TO=1] ~ ~ ~
/execute @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=2,score_ACV_COLOR=2] ~ ~ ~ teleport @e[score_ACV_TP_TO_min=2,score_ACV_TP_TO=2] ~ ~ ~
/execute @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=3,score_ACV_COLOR=3] ~ ~ ~ teleport @e[score_ACV_TP_TO_min=3,score_ACV_TP_TO=3] ~ ~ ~

// Rotate the entity
/tp @e[score_ACV_TP_RD_min=0,score_ACV_TP_RD=0] ~ ~ ~ ~180 ~
/tp @e[score_ACV_TP_RD_min=1,score_ACV_TP_RD=1] ~ ~ ~ ~90 ~
//tp @e[score_ACV_TP_RD_min=2,score_ACV_TP_RD=2] ~ ~ ~ ~0 ~
/tp @e[score_ACV_TP_RD_min=3,score_ACV_TP_RD=3] ~ ~ ~ ~-90 ~

// when IR is up or down, but OR is on a wall (or it is a lightbridge)
/tp @e[score_ACV_TP_IR=-1,score_ACV_TP_OR_min=0,score_ACV_TP_OR=0] ~ ~ ~ 0 ~
/tp @e[score_ACV_TP_IR=-1,score_ACV_TP_OR_min=1,score_ACV_TP_OR=1] ~ ~ ~ 90 ~
/tp @e[score_ACV_TP_IR=-1,score_ACV_TP_OR_min=2,score_ACV_TP_OR=2] ~ ~ ~ 180 ~
/tp @e[score_ACV_TP_IR=-1,score_ACV_TP_OR_min=3,score_ACV_TP_OR=3] ~ ~ ~ -90 ~

// Lightbridges must look also look away from horizontal portals
/tp @e[tag=ACV_RotationProof,score_ACV_TP_OR_min=-1,score_ACV_TP_OR=-1] ~ ~ ~ ~ -90
/tp @e[tag=ACV_RotationProof,score_ACV_TP_OR_min=-2,score_ACV_TP_OR=-2] ~ ~ ~ ~ 90
}


text: sort([3,2,4,5],output = permutation);


text:
######################################################################

`is_element/cubes` := (k::posint) -> (A::set) -> proc(f)
 local n,i,j;
 global reason;

 if not type(f,table) then
  reason := [convert(procname,string),"f is not a table",f];
  return false;
 fi;

 if map(op,{indices(f)}) <> A then
  reason := [convert(procname,string),"f is not indexed by A",f,A];
  return false;
 fi;

 n := nops(A);
 for i from 1 to n do
  if not `is_element/single_cubes`(k)(f[A[i]]) then
   reason := [convert(procname,string),"f(a) is not a single cube",A[i],f(A[i]),reason];
   return false;
  fi;

  for j from 1 to i-1 do
   if `overlap/single_cubes`(k)(f[A[j]],f[A[i]]) then
    reason := [convert(procname,string),"subcubes f(a) and f(b) overlap",A[i],A[j],f(A[i]),f(A[j])];
    return false;
   fi;
  od;
 od;

 return true;
end;

`is_equal/cubes` := (k::posint) -> (A::set) -> proc(f,g)
 local a;
 global reason;

 for a in A do
  if not `is_equal/single_cubes`(k)(f[a],g[a]) then
   reason := [convert(procname,string),"f[a] <> g[a]",a,f[a],g[a]];
   return false;
  fi;
 od;

 return true;
end;

`is_leq/cubes` := NULL;

`random_element/cubes` := (N::posint) -> (A::set) -> proc()
 local d,f,a;

 d := 5;

 while true do
  d := d+1;
  f := table();
  for a in A do
   f[a] := `random_element/single_cubes`(N)(d);
  od;

  if `is_element/cubes`(N)(A)(f) then
   return(eval(f));
  fi;

  if d > 1000 then return FAIL; fi;
 od;
end;

`list_elements/cubes` := NULL;
`count_elements/cubes` := NULL;

`draw/cubes` := (A::set) -> proc(f)
 local P,a;

 P := rectangle([0,0],[1,1],colour=blue,style=line);
 for a in A do
  P := P,rectangle(f[a][1],f[a][2],colour=red);
 od;
 display(P,scaling=constrained,axes=none);
end;

`centres/cubes` := (k::posint) -> (A::set) -> proc(f)
 local g,a;

 g := table();
 for a in A do
  g[a] := `centre/single_cubes`(k)(f[a]);
 od;

 return eval(g);
end;


text:
# This code sets up a toy model of the nonstandard product structure
# on the HKR character ring for V x V, where V is the groupoid of
# finite-dimensional F-vector spaces.

# We fix a number d > 0 (which should really be infinity).
# We consider the set A1 of tables indexed by natural numbers u < d,
# and the set A2 indexed by pairs (u,v) with u+v < d.  These sets
# have products called `star/A1` and `star/A2`.  The first is given by
# the following rule: `star/A1`(f,g)[u] is the sum of all terms
# c1[u0,u1] * f[u0] * g[u1], where u0+u1 = u, and c1[u0,u1] is the
# number of splittings of F^u as a direct sum of subspaces of
# dimensions u0 and u1.
#
# Similarly, `star/A2`(f,g)[u,v] is a sum of terms
# c2[u0,v0,u1,v1] * f[u0,v0] * g[u1,v1] where
# (u0,v0) + (u1,v1) = (u,v) and
# c2[u0,v0,u1,v1] is c1[u0,u1] times c1[v0,v1] times
# q^(u0*v1+2*u1*v0).  This final factor is inserted to ensure that
# the coproduct map psi : A1 -> A2 (given by psi(f)[u,v]=f[u+v])
# is a ring map.

d := 10;

grassmann_count := (q) -> (n,d) ->
 mul(q^(n-i)-1,i=0..d-1)/mul(q^(d-i)-1,i=0..d-1);

splitting_count := (q) -> (n,m) ->
 grassmann_count(q)(n+m,n) * q^(n*m);

######################################################################

`is_element/A1` := proc(f)
 local u,L;

 if not(type(f,table)) then
  return false;
 fi;

 L := [seq([u],u=0..d-1)];

 if [indices(f)] <> L then
  return false;
 fi;

 return true;
end:

`is_equal/A1` := proc(f,g)
 local u;

 for u from 0 to d-1 do
  if simplify(f[u]-g[u]) <> 0 then
   return false;
  fi;
 od;
 return true;
end:

`to_list/A1` := (f) -> [seq(f[u],u=0..d-1)]:
`from_symbol/A1` := (f) -> table([seq(u=f[u],u=0..d-1)]);

`zero/A1`  := table([seq(u=0,u=0..d-1)]):
`one/A1`   := table([seq(u=1,u=0..d-1)]):
`delta/A1` := table([seq(u=`if`(u=0,1,0),u=0..d-1)]):

`plus/A1` := proc(f,g)
 table([seq(u = f[u]+g[u],u=0..d-1)]);
end:

`minus/A1` := proc(f,g)
 table([seq(u = f[u]-g[u],u=0..d-1)]);
end:

`dot/A1` := proc(f,g)
 table([seq(u = f[u]*g[u],u=0..d-1)]);
end:

`scalar_times/A1` := proc(c,f)
 table([seq(u = c*f[u],u=0..d-1)]);
end:

`star/A1` := proc(f,g)
 local h,u,u0,u1;
 h := table();
 for u from 0 to d-1 do
  h[u] := 0;
 od;

 for u0 from 0 to d-1 do
  for u1 from 0 to d-1-u0 do
   h[u0+u1] := h[u0+u1] + splitting_count(q)(u0,u1)*f[u0]*g[u1];
  od;
 od:

 return eval(h);
end:

######################################################################

`is_element/A2` := proc(f)
 local u,L;

 if not(type(f,table)) then
  return false;
 fi;

 L := [seq(seq([u,v],v=0..d-1-u),u=0..d-1)];

 if [indices(f)] <> L then
  return false;
 fi;

 return true;
end:

`is_equal/A2` := proc(f,g)
 local u,v;

 for u from 0 to d-1 do
  for v from 0 to d-1-u do
   if simplify(f[u,v]-g[u,v]) <> 0 then
    return false;
   fi;
  od;
 od;
 return true;
end:

`to_list/A2` := (f) -> [seq(seq(f[u,v],v=0..d-1-u),u=0..d-1)]:
`from_symbol/A2` := (f) -> table([seq(seq((u,v)=f[u,v],v=0..d-1-u),u=0..d-1)]);

`zero/A2`  := table([seq(seq((u,v)=0,v=0..d-1-u),u=0..d-1)]):
`one/A2`   := table([seq(seq((u,v)=1,v=0..d-1-u),u=0..d-1)]):
`delta/A2` := table([seq(seq((u,v)=`if`(u=0 and v = 0,1,0),v=0..d-1-u),u=0..d-1)]):

`plus/A2` := proc(f,g)
 table([seq(seq((u,v) = f[u,v]+g[u,v],v=0..d-1-u),u=0..d-1)]);
end:

`minus/A2` := proc(f,g)
 table([seq(seq((u,v) = f[u,v]-g[u,v],v=0..d-1-u),u=0..d-1)]);
end:

`dot/A2` := proc(f,g)
 table([seq(seq((u,v) = f[u,v]*g[u,v],v=0..d-1-u),u=0..d-1)]);
end:

`scalar_times/A2` := proc(c,f)
 table([seq(seq((u,v) = c*f[u,v],v=0..d-1-u),u=0..d-1)]);
end:

`mu/A2` := (u0,v0,u1,v1) -> q^(u0*v1+2*u1*v0);

`star/A2` := proc(f,g)
 local h,u,v,u0,u1,v0,v1;
 h := table();

 for u from 0 to d-1 do
  for v from 0 to d-1-u do
   h[u,v] := 0;
  od;
 od:

 for u0 from 0 to d-1 do
  for u1 from 0 to d-1-u0 do
   for v0 from 0 to d-1-u0-u1 do
    for v1 from 0 to d-1-u0-u1-v0 do
     h[u0+u1,v0+v1] :=
      h[u0+u1,v0+v1] +
       splitting_count(q)(u0,u1) *
       splitting_count(q)(v0,v1) *
       `mu/A2`(u0,v0,u1,v1) *
       f[u0,v0] * g[u1,v1];
    od;
   od;
  od;
 od;

 return eval(h);
end:

######################################################################

`is_element/A3` := proc(f)
 local u,L;

 if not(type(f,table)) then
  return false;
 fi;

 L := [seq(seq(seq([u,v,w],w=0..d-1-u-v),v=0..d-1-u),u=0..d-1)];

 if [indices(f)] <> L then
  return false;
 fi;

 return true;
end:

`is_equal/A3` := proc(f,g)
 local u,v,w;

 for u from 0 to d-1 do
  for v from 0 to d-1-u do
   for w from 0 to d-1-u-v do
    if simplify(f[u,v,w]-g[u,v,w]) <> 0 then
     return false;
    fi;
   od;
  od;
 od;
 return true;
end:

`to_list/A3` := (f) -> [seq(seq(seq(f[u,v,w],w=0..d-1-u-v),v=0..d-1-u),u=0..d-1)]:
`from_symbol/A3` := (f) -> table([seq(seq(seq((u,v,w)=f[u,v,w],w=0..d-1-u-v),v=0..d-1-u),u=0..d-1)]);

`zero/A3`  := table([seq(seq(seq((u,v,w)=0,w=0..d-1-u-v),v=0..d-1-u),u=0..d-1)]):
`one/A3`   := table([seq(seq(seq((u,v,w)=1,w=0..d-1-u-v),v=0..d-1-u),u=0..d-1)]):
`delta/A3` := table([seq(seq(seq((u,v,w)=`if`([u,v,w]=[0,0,0],1,0),w=0..d-1-u-v),v=0..d-1-u),u=0..d-1)]):

`plus/A3` := proc(f,g)
 table([seq(seq(seq((u,v,w) = f[u,v,w]+g[u,v,w],w=0..d-1-u-v),v=0..d-1-u),u=0..d-1)]);
end:

`minus/A3` := proc(f,g)
 table([seq(seq(seq((u,v,w) = f[u,v,w]-g[u,v,w],w=0..d-1-u-v),v=0..d-1-u),u=0..d-1)]);
end:

`dot/A3` := proc(f,g)
 table([seq(seq(seq((u,v,w) = f[u,v,w]*g[u,v,w],w=0..d-1-u-v),v=0..d-1-u),u=0..d-1)]);
end:

`scalar_times/A3` := proc(c,f)
 table([seq(seq(seq((u,v,w) = c*f[u,v,w],w=0..d-1-u-v),v=0..d-1-u),u=0..d-1)]);
end:

`mu/A3` := (u0,v0,w0,u1,v1,w1) -> q^(u0*v1+u0*w1+v0*w1+2*u1*v0+2*u1*w0+2*v1*w0);

`star/A3` := proc(f,g)
 local h,u,v,w,u0,u1,v0,v1,w0,w1;
 h := table();

 for u from 0 to d-1 do
  for v from 0 to d-1-u do
   for w from 0 to d-1-u-v do
    h[u,v,w] := 0;
   od;
  od;
 od:

 for u0 from 0 to d-1 do
  for u1 from 0 to d-1-u0 do
   for v0 from 0 to d-1-u0-u1 do
    for v1 from 0 to d-1-u0-u1-v0 do
     for w0 from 0 to d-1-u0-u1-v0-v1 do
      for w1 from 0 to d-1-u0-u1-v0-v1-w0 do
       h[u0+u1,v0+v1,w0+w1] :=
	h[u0+u1,v0+v1,w0+w1] +
	 splitting_count(q)(u0,u1) *
	 splitting_count(q)(v0,v1) *
	 splitting_count(q)(w0,w1) *
	 `mu/A3`(u0,v0,w0,u1,v1,w1) *
	 f[u0,v0,w0] * g[u1,v1,w1];
      od;
     od;
    od;
   od;
  od;
 od;

 return eval(h);
end:

######################################################################

psi := proc(f)
 local h,u,v;

 h := table();

 for u from 0 to d-1 do
  for v from 0 to d-1-u do
   h[u,v] := f[u+v];
  od;
 od:

 return eval(h);
end:

psi_L := proc(f)
 local h,u,v,w;

 h := table();

 for u from 0 to d-1 do
  for v from 0 to d-1-u do
   for w from 0 to d-1-u-v do
    h[u,v,w] := f[u+v,w];
   od;
  od;
 od:

 return eval(h);
end:

psi_R := proc(f)
 local h,u,v,w;

 h := table();

 for u from 0 to d-1 do
  for v from 0 to d-1-u do
   for w from 0 to d-1-u-v do
    h[u,v,w] := f[u,v+w];
   od;
  od;
 od:

 return eval(h);
end:

tensor11 := proc(f,g)
 table([seq(seq((u,v)=f[u]*g[v],v=0..d-1-u),u=0..d-1)]);
end:

twisted_tensor11 := proc(f,g)
 table([seq(seq((u,v)=q^(u*v)*f[u]*g[v],v=0..d-1-u),u=0..d-1)]);
end:

eta_L := (f) -> tensor11(`delta/A1`,f):

eta_R := (f) -> tensor11(f,`delta/A1`):

epsilon_L := proc(f)
 local h,u;

 h := table;

 for u from 0 to d-1 do
  h[u] := f[0,u];
 od;

 return eval(h):
end:

epsilon_R := proc(f)
 local h,u;

 h := table;

 for u from 0 to d-1 do
  h[u] := f[u,0];
 od;

 return eval(h):
end:

######################################################################

check_A1 := proc()
 local f,g,h,ff,gg,hh;

 ff := `from_symbol/A1`(f);
 gg := `from_symbol/A1`(g);
 hh := `from_symbol/A1`(h);

 printf("%a()\n",procname);

 _ASSERT(
  `is_equal/A1`(`star/A1`(`delta/A1`,ff),ff),"Left unit law"
 );

 _ASSERT(
  `is_equal/A1`(`star/A1`(ff,`delta/A1`),ff),"Right unit law"
 );

 _ASSERT(
  `is_equal/A1`(`star/A1`(ff,`star/A1`(gg,hh)),
                `star/A1`(`star/A1`(ff,gg),hh)),
  "Associativity law"
 );
end:

check_A2 := proc()
 local f,g,h,ff,gg,hh;

 ff := `from_symbol/A2`(f);
 gg := `from_symbol/A2`(g);
 hh := `from_symbol/A2`(h);

 printf("%a()\n",procname);

 _ASSERT(
  `is_equal/A2`(`star/A2`(`delta/A2`,ff),ff),"Left unit law"
 );

 _ASSERT(
  `is_equal/A2`(`star/A2`(ff,`delta/A2`),ff),"Right unit law"
 );

 _ASSERT(
  `is_equal/A2`(`star/A2`(ff,`star/A2`(gg,hh)),
                `star/A2`(`star/A2`(ff,gg),hh)),
  "Associativity law"
 );
end:

check_A3 := proc()
 local f,g,h,ff,gg,hh;

 ff := `from_symbol/A3`(f);
 gg := `from_symbol/A3`(g);
 hh := `from_symbol/A3`(h);

 printf("%a()\n",procname);

 _ASSERT(
  `is_equal/A3`(`star/A3`(`delta/A3`,ff),ff),"Left unit law"
 );

 _ASSERT(
  `is_equal/A3`(`star/A3`(ff,`delta/A3`),ff),"Right unit law"
 );

 _ASSERT(
  `is_equal/A3`(`star/A3`(ff,`star/A3`(gg,hh)),
                `star/A3`(`star/A3`(ff,gg),hh)),
  "Associativity law"
 );
end:

check_A_maps := proc()
 local f1,g1,h1,ff1,gg1,hh1,f2,g2,h2,ff2,gg2,hh2;

 ff1 := `from_symbol/A1`(f1);
 gg1 := `from_symbol/A1`(g1);
 hh1 := `from_symbol/A1`(h1);

 ff2 := `from_symbol/A2`(f2);
 gg2 := `from_symbol/A2`(g2);
 hh2 := `from_symbol/A2`(h2);

 _ASSERT(
  `is_equal/A2`(eta_L(`delta/A1`),`delta/A2`) and
  `is_equal/A2`(eta_L(`star/A1`(ff1,gg1)),`star/A2`(eta_L(ff1),eta_L(gg1))),
  "eta_L is a ring map"
 );

 _ASSERT(
  `is_equal/A2`(eta_R(`delta/A1`),`delta/A2`) and
  `is_equal/A2`(eta_R(`star/A1`(ff1,gg1)),`star/A2`(eta_R(ff1),eta_R(gg1))),
  "eta_R is a ring map"
 );

 _ASSERT(
  `is_equal/A2`(psi(`delta/A1`),`delta/A2`) and
  `is_equal/A2`(psi(`star/A1`(ff1,gg1)),`star/A2`(psi(ff1),psi(gg1))),
  "psi is a ring map"
 );

 _ASSERT(
  `is_equal/A1`(epsilon_L(`delta/A2`),`delta/A1`) and
  `is_equal/A1`(epsilon_L(`star/A2`(ff2,gg2)),`star/A1`(epsilon_L(ff2),epsilon_L(gg2))),
  "epsilon_L is a ring map"
 );

 _ASSERT(
  `is_equal/A1`(epsilon_R(`delta/A2`),`delta/A1`) and
  `is_equal/A1`(epsilon_R(`star/A2`(ff2,gg2)),`star/A1`(epsilon_R(ff2),epsilon_R(gg2))),
  "epsilon_R is a ring map"
 );

 _ASSERT(
  `is_equal/A1`(epsilon_L(psi(ff1)),ff1),
  "epsilon_L o psi = 1"
 );

 _ASSERT(
  `is_equal/A1`(epsilon_R(psi(ff1)),ff1),
  "epsilon_R o psi = 1"
 );

 _ASSERT(
  `is_equal/A2`(`star/A2`(eta_R(ff1),eta_L(gg1)),twisted_tensor11(ff1,gg1)),
  "eta_R * eta_L"
 );
end;


text:
make_fgl_Lazard := proc(d)
 global a,m,ma,am;
 local T,u,x,y,q,i,k,Fm,sol;

 T := table();
 T["degree"] := d;
 Order := d+2;

 T["log_m"] := unapply(x + add(m[i]*x^(i+1),i=1..d-1),x);
 T["exp_m"] :=
   unapply(
   convert(
    solve(x = series(T["log_m"](y) + sin(y)^(d+1),y=0,d+1),y),
    polynom,x
   ),
   x
  );

 T["sum_m"] :=
  unapply(
   subs(u=1,expand(convert(series(T["exp_m"](T["log_m"](u*x)+T["log_m"](u*y)),u=0,d+1),polynom,u))),
   x,y
  ):

 Fm := T["sum_m"](x,y);
 am[1] := m[1];
 for k from 2 to d-1 do
  q := igcd_alt([seq(binomial(k+1,i),i=1..k)])[2];
  am[k] := -add(q[i]*coeff(coeff(Fm,x,i),y,k+1-i),i=1..k);
 od:
 sol := solve({seq(a[i]=am[i],i=1..d-1)},{seq(m[i],i=1..d-1)}):
 for i from 1 to d-1 do ma[i] := subs(sol,m[i]); od:

 T["log_a"] := unapply(collect(expand(eval(subs(m = ma,T["log_m"](x)))),x),x);
 T["exp_a"] := unapply(collect(expand(eval(subs(m = ma,T["exp_m"](x)))),x),x);
 T["sum_a"] := unapply(collect(expand(eval(subs(m = ma,T["sum_m"](x,y)))),{x,y}),x,y);

 return eval(T);
end:


text:
// Setup
include "ACV_setupScoreboards.mpl"
include "ACV_deleteScoreboards.mpl"

// Main
include "ACV_portals.mpl"

// Remote
include "ACV_remote_deletePortals.mpl"
include "ACV_remote_deleteAllPortals.mpl"
include "ACV_remote_createPortals.mpl"

// Extra
include "ACV_confirmBlockOption.mpl"
include "ACV_teleportation.mpl"
include "ACV_smoothTeleportation.mpl"
include "ACV_lightBridges.mpl"
include "ACV_antiBridges.mpl"
include "ACV_cubes.mpl"
include "ACV_lasers.mpl"
include "ACV_megs.mpl"
include "ACV_gels.mpl"

project ACV {
  orientation "zyx"
}

install {
/gamerule commandBlockOutput false
/gamerule logAdminCommands false
/tellraw @a {"text":"[ApertureCraft Vanilla] ","color":"red","extra":[{"text":"Starting Installation","color":"yellow"}]}

/tellraw @a {"text":"[ApertureCraft Vanilla] ","color":"gold","extra":[{"text":"Setting up Scoreboards","color":"yellow"}]}
ACV_setupScoreboards()
/tellraw @a {"text":"[ApertureCraft Vanilla] ","color":"gold","extra":[{"text":"Enabling UI","color":"yellow"}]}
/fill ${origin + (9 1 -7)} ${origin + (15 3 -6)} lever facing=north replace barrier
/blockdata ${origin + (7 1 -16)} {Command:"/setblock ~-7 ~ ~16 redstone_block"}
/setblock ${origin + (7 2 -14)} wall_sign 3 replace {Text1:"{\"text\":\"\"}",Text2:"{\"text\":\"[Uninstall]\",\"color\":\"dark_green\",\"bold\":\"true\"}",Text3:"{\"text\":\"\"}",Text4:"{\"text\":\"\"}"}
/fill ${origin + (6 0 -15)} ${origin + (9 2 -15)} stained_hardened_clay 13 replace stained_hardened_clay 14

/tellraw @a {"text":"[ApertureCraft Vanilla] ","color":"gold","extra":[{"text":"Setting up Block-Option","color":"yellow"}]}
ACV_confirmBlockOption()
}

uninstall {
/tellraw @a {"text":"[ApertureCraft Vanilla] ","color":"red","extra":[{"text":"Starting Uninstallation","color":"yellow"}]}

/tellraw @a {"text":"[ApertureCraft Vanilla] ","color":"gold","extra":[{"text":"Deleting Portals","color":"yellow"}]}
ACV_remote_deleteAllPortals()

/tellraw @a {"text":"[ApertureCraft Vanilla] ","color":"gold","extra":[{"text":"Disabling UI","color":"yellow"}]}
/fill ${origin + (9 1 -7)} ${origin + (15 3 -6)} barrier 0 replace lever
/blockdata ${origin + (7 1 -16)} {Command:"/setblock ~-7 ~-1 ~16 redstone_block"}
/setblock ${origin + (7 2 -14)} wall_sign 3 replace {Text1:"{\"text\":\"\"}",Text2:"{\"text\":\"[Install]\",\"color\":\"dark_red\",\"bold\":\"true\"}",Text3:"{\"text\":\"\"}",Text4:"{\"text\":\"\"}"}
/fill ${origin + (6 0 -15)} ${origin + (9 2 -15)} stained_hardened_clay 14 replace stained_hardened_clay 13

///tellraw @a {"text":"[ApertureCraft Vanilla] ","color":"gold","extra":[{"text":"Interrupt all Interface calls","color":"yellow"}]}
///execute @e[tag=ACV_remote_interface] ~ ~ ~ /setblock ~ ~ ~ stone
///execute @e[tag=ACV_internal_interface] ~ ~ ~ /setblock ~ ~ ~ stone

/tellraw @a {"text":"[ApertureCraft Vanilla] ","color":"gold","extra":[{"text":"Deleting Scoreboards","color":"yellow"}]}
ACV_deleteScoreboards()

// 3 tick delay
/setblock ${this+1} redstone_block
skip
impulse: /setblock ${this-1} stone
/setblock ${this+1} redstone_block
skip
impulse: /setblock ${this-1} stone
/setblock ${this+1} redstone_block
skip
impulse: /setblock ${this-1} stone
// 3 tick delay end

/tellraw @a {"text":"[ApertureCraft Vanilla] ","color":"gold","extra":[{"text":"Deleting Interfaces","color":"yellow"}]}
}

max_stars_repo_path: src/main/mpl/ApertureCraft Vanilla.mpl

max_issues_repo_path: src/main/mpl/ApertureCraft Vanilla.mpl

max_forks_repo_path: src/main/mpl/ApertureCraft Vanilla.mpl


text:
# Teach Maple (through depends and eval) about our new binding forms.
# Integrand and LO bind from 1st arg to 2nd arg.

`depends/Integrand` := proc(v, e, x) depends(e, x minus {v}) end proc:
`depends/LO`        := proc(v, e, x) depends(e, x minus {v}) end proc:

`eval/Integrand` := proc(e, eqs)
  local v, ee;
  v, ee := op(e);
  eval(op(0,e), eqs)(BindingTools:-generic_evalat(v, ee, eqs))
end proc:

`eval/LO` := proc(e, eqs)
  local v, ee;
  v, ee := op(e);
  eval(op(0,e), eqs)(BindingTools:-generic_evalat(v, ee, eqs))
end proc:

#############################################################################

NewSLO := module ()
  option package;
  local
        t_sum, t_product,
        mysolve, Shiftop, Diffop, Recognized,
        bind, weight,
        reduce_IntsSums, reduce_Integrals,
        elim_intsum, int_assuming, sum_assuming,
        banish, isBound_IntsSums, isBound_IntSum,
        mk_sym, mk_ary, mk_idx, innermostIntSum, ChangeVarInt, SimplifyKB_,
        ModuleLoad;
  export
     # These first few are smart constructors (for themselves):
         integrate, applyintegrand,
     # while these are "proper functions"
         RoundTrip, Simplify, SimplifyKB,  apply_LO,
         TestSimplify, TestHakaru, TestDisint, Efficient, TestEfficient,
         Concrete,
         toLO, fromLO, improve, reduce,
         density, bounds, unweight,

     # Like simplify_assuming, but also applies `hack_Beta` and `eval_factor`
     # which helps some things simplify.
         simplify_factor_assuming,

         ReparamDetermined, determined, reparam, disint ;
  # these names are not assigned (and should not be).  But they are
  # used as global names, so document that here.
  global LO, Integrand, SumIE, ProductIE;
  uses Hakaru, KB, Loop, Partition, Domain, disint;

  t_sum     := 'specfunc({sum    ,Sum    })';
  t_product := 'specfunc({product,Product})';

$include "NewSLO/Interface.mpl"
$include "NewSLO/To.mpl"
$include "NewSLO/From.mpl"
$include "NewSLO/Banish.mpl"
$include "NewSLO/Improve.mpl"
$include "NewSLO/Disint.mpl"
$include "NewSLO/Reparam.mpl"
$include "NewSLO/Factor.mpl"

  # This is for backwards compatibility; do not use them directly,
  # instead use the symbols exported from inner modules of NewSLO
  `integrate` := toLO:-integrate;
  `unweight`  := fromLO:-unweight;

  apply_LO := proc(e::specfunc(LO),f,$) unapply(op(2,e),op(1,e))(f) end proc;

# An integrand h is either an Integrand (our own binding construct for a
# measurable function to be integrated) or something that can be applied
# (probably proc, which should be applied immediately, or a generated symbol).
  applyintegrand := proc(h, x, $)
    local var, body, dummy;
    if h :: 'Integrand(name, anything)' then
      var, body := op(h);
      if x :: {boolean, specfunc({And,Or,Not})} then
        body := eval(body, var=dummy);
        var  := dummy;
        body := subs((var=true)=var, body);
      end if;
      eval(body, var=x)
    elif h :: appliable then
      h(x)
    else
      'procname(_passed)'
    end if
  end proc;

  mysolve := proc(constraints)
    # This wrapper around "solve" works around the problem that Maple sometimes
    # thinks there is no solution to a set of constraints because it doesn't
    # recognize the solution to each constraint is the same.  For example--
    # This fails     : solve({c*2^(-1/2-alpha) = sqrt(2)/2, c*4^(-alpha) = 2^(-alpha)}, {c}) assuming alpha>0;
    # This also fails: solve(simplify({c*2^(-1/2-alpha) = sqrt(2)/2, c*4^(-alpha) = 2^(-alpha)}), {c}) assuming alpha>0;
    # But this works : map(solve, {c*2^(-1/2-alpha) = sqrt(2)/2, c*4^(-alpha) = 2^(-alpha)}, {c}) assuming alpha>0;
    # And the difference of the two solutions returned simplifies to zero.

    local result;
    if nops(constraints) = 0 then return NULL end if;
    result := solve(constraints, _rest);
    if result <> NULL or not (constraints :: {set,list}) then
      return result
    end if;
    result := mysolve(subsop(1=NULL,constraints), _rest);
    if result <> NULL
       and op(1,constraints) :: 'anything=anything'
       and simplify(eval(op([1,1],constraints) - op([1,2],constraints),
                         result)) <> 0 then
      return NULL
    end if;
    result
  end proc;

  ###
  # smart constructors for our language

  bind := proc(m, x, n, $)
    if n = 'Ret'(x) then
      m # monad law: right identity
    elif m :: 'Ret(anything)' then
      eval(n, x = op(1,m)) # monad law: left identity
    else
      'Bind(_passed)'
    end if;
  end proc;

  weight := proc(p, m, $)
    #Trying to make the below into an ASSERT statement results in a kernel
    #crash, even though the condition is caught, the message printed, and
    #the error generated.
    if kernelopts(assertlevel) > 0 then
         if not p::
              (complexcons &implies ((numeric &implies {0,1}) &under csgn))
         then
              userinfo(
                   1, SLO,
                   proc()
                   uses LL= ListTools:-LengthSplit;
                   local
                        S:= Matrix([LL(debugopts(callstack)[2..], 3)]),
                        rts:= interface(rtablesize= infinity)
                   ;
                        print(S);
                        interface(rtablesize= rts);
                        NULL
                   end proc()
              );
              error "Negative weight %1 not allowed", p
         end if
    end if;
    if p = 1 then
      m
    elif p = 0 then
      Msum()
    elif m :: 'Weight(anything, anything)' then
      weight(p * op(1,m), op(2,m))
    else
      'Weight(_passed)'
    end if;
  end proc;

  density[Lebesgue] := proc(lo,hi,$) proc(x,$) 1 end proc end proc;
  density[Uniform] := proc(a,b,$) proc(x,$)
    1/(b-a)
  end proc end proc;
  density[Gaussian] := proc(mu, sigma, $) proc(x,$)
    1/sigma/sqrt(2)/sqrt(Pi)*exp(-(x-mu)^2/2/sigma^2)
  end proc end proc;
  density[Cauchy] := proc(loc, scale, $) proc(x,$)
    1/Pi/scale/(1+((x-loc)/scale)^2)
  end proc end proc;
  density[StudentT] := proc(nu, loc, scale, $) proc(x,$)
    GAMMA((nu+1)/2) / GAMMA(nu/2) / sqrt(Pi*nu) / scale
    * (1 + ((x-loc)/scale)^2/nu)^(-(nu+1)/2)
  end proc end proc;
  density[BetaD] := proc(a, b, $) proc(x,$)
    x^(a-1)*(1-x)^(b-1)/Beta(a,b)
  end proc end proc;
  # Hakaru uses the alternate definition of gamma, so the args are backwards
  density[GammaD] := proc(shape, scale, $) proc(x,$)
    x^(shape-1)/scale^shape*exp(-x/scale)/GAMMA(shape);
  end proc end proc;
  density[InverseGammaD]:= proc(shape, scale, $)
       proc(x,$) scale^shape/GAMMA(shape)/x^(shape+1)/exp(scale/x) end proc
  end proc;
  density[Counting] := proc(lo, hi, $) proc(k,$)
    1
  end proc end proc;
  density[Categorical] := proc(a, $) proc(k,$)
    idx(a,k)
  end proc end proc;
  density[Binomial]:= proc(n,p,$)
       proc(k,$) p^k*(1-p)^(n-k)*GAMMA(n+1)/GAMMA(k+1)/GAMMA(n-k+1) end proc
  end proc;
  density[NegativeBinomial] := proc(r, p, $) proc(k,$)
    p^k * (1-p)^r * GAMMA(r+k) / GAMMA(k+1) / GAMMA(r)
  end proc end proc;
  density[PoissonD] := proc(lambda, $) proc(k,$)
    lambda^k/exp(lambda)/k!
  end proc end proc;
  density[ChiSquared] := proc(k, $) proc(x,$)
    x^((1/2)*k-1)*exp(-(1/2)*x)/(2^((1/2)*k)*GAMMA((1/2)*k))
  end proc end proc:

  bounds[Lebesgue] := `..`;
  bounds[Uniform] := proc(a, b, $) a .. b end proc;
  bounds[Gaussian] := proc(mu, sigma, $) -infinity .. infinity end proc;
  bounds[Cauchy] := proc(loc, scale, $) -infinity .. infinity end proc;
  bounds[StudentT] := proc(nu, loc, scale, $) -infinity .. infinity end proc;
  bounds[BetaD] := proc(a, b, $) 0 .. 1 end proc;
  bounds[GammaD] := proc(shape, scale, $) 0 .. infinity end proc;
  bounds[InverseGammaD]:= proc(shape, scale, $) 0..infinity end proc;
  bounds[Counting] := proc(lo, hi, $) lo..hi-1 end proc;
  bounds[Categorical] := proc(a, $) 0 .. size(a)-1 end proc;
  bounds[Binomial]:= proc(n,p,$) 0..n end proc;
  bounds[NegativeBinomial] := proc(r, p, $) 0 .. infinity end proc;
  bounds[PoissonD] := proc(lambda, $) 0 .. infinity end proc;
  bounds[ChiSquared] := proc(k, $) 0 .. infinity end proc;

  mk_sym := proc(var :: name, h, $)
    local x;
    x := var;
    if h :: 'Integrand(name, anything)' then
      x := op(1,h);
    elif h :: 'procedure' then
      x := op(1, [op(1,h), x]);
    end if;
    gensym(x)
  end proc;

  mk_ary := proc(e, loops :: list(name = range), $)
    foldl((res, i) -> ary(op([2,2],i) - op([2,1],i) + 1,
                          op(1,i),
                          eval(res, op(1,i) = op(1,i) + op([2,1],i))),
          e, op(loops));
  end proc;

  mk_idx := proc(e, loops :: list(name = range), $)
    foldr((i, res) -> idx(res, op(1,i) - op([2,1],i)),
          e, op(loops));
  end proc;

  ModuleLoad := proc($)
    local prev, ex;
    ModuleUnload();
    Hakaru; # Make sure the Hakaru module is loaded for the types t_type and t_Hakaru.
    KB;     # Make sure the KB module is loaded, for the type t_kb
    prev := kernelopts(opaquemodules=false);
    try
      PiecewiseTools:-InertFunctions := PiecewiseTools:-InertFunctions
        union '{# Do not lift piecewise over a binder
                Integrand,LO,lam,Branch,Bind,ary,Plate,
                forall,Ints,Sums,ints,sums,`..`}';
    catch:
         userinfo(
              1, NewSLO,
              "Redefinition of PiecewiseTools:-InertFunctions failed.",
              StringTools:-FormatMessage(lastexception[2..-1])
         )
    finally
      kernelopts(opaquemodules=prev);
    end try;

    for ex in [ [`Int`,isBound_IntSum], [`Sum`,isBound_IntSum]
              , [`Ints`,isBound_IntsSums], [`Sums`,isBound_IntsSums] ] do
      Domain:-Set_ExtBound(op(1,ex), op(2,ex)(op(1,ex)));
    end do;

    for ex in [ [`PARTITION`, Partition:-isShape] , [ 'piecewise', Partition:-isShape] ] do
      Domain:-Set_ExtShape(op(1,ex), op(2,ex)(op(1,ex)));
    end do;

  end proc; #ModuleLoad
  ModuleLoad():
end module; # NewSLO


text:
# Kinematic constraints for Cloos QRC 350
# Description
# The robot has a coupling at the hand axes. This worksheet is necessary to fit this into the dynamics toolbox.
# The constraints can also be declared solely in the robot_env. This however leads to some functions not working as expected
# Sources
# [Diekmeyer2018_S678] Diekmeyer, Jonas: Identifikation der inversen Dynamik eines seriellen Roboters im geschlossenen Regelkreis; Studienarbeit; Leibniz Universität Hannover
#
# Author
# Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2020-06
# (C) Institut für Mechatronische Systeme, Leibniz Universität Hannover
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
codegen_act := true:
codegen_opt := 2: # Hoher Optimierungsgrad.
;
read "../helper/proc_MatlabExport":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
# Variable mit Winkeln der Zwangsbedingungen nur in Abhängigkeit der verallgemeinerten Koordinaten
kintmp_qs := Matrix(RowDimension(kintmp_s),1):
kintmp_qt := Matrix(RowDimension(kintmp_s),1):
# Ersetzungsausdrücke definieren.
# Variable zum speichern des Sinus und Cosinus der Winkel. Für dieses System ist das eigentlich nicht notwendig. Erstelle Variable, da sie von den anderen Skripten erwartet wird
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
# Kinematic Constraints
# First axis is defined in negative direction
rho1_qt := -qJ_t(1,1):
# Hand axes 5 and 6 have a coupling via differential gear.
# [Diekmeyer2018_S678], p. 24
rho6_qt := qJ_t(6,1) - kDG*qJ_t(5,1):
kintmp_qt(1,1) := rho1_qt:
kintmp_qt(2,1) := rho6_qt:
# Nachverarbeitung
# Umrechnung in Substitutionsvariablen
kintmp_qs := convert_t_s(kintmp_qt):
# Exportiere Code für folgende Skripte
# Speichere Maple-Ausdruck (Eingabe-Format und internes Format)
save kintmp_qs, kintmp_qt, kin_constraints_exist, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert", robot_name):
save kintmp_qs, kintmp_qt, kin_constraints_exist, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_qs )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(Transpose(kc_symbols), sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m", robot_name), 2):