1 /* MetaEthical.AI
 2    ==============
 3    Read the Introduction at:
 4    https://medium.com/@june_ku/formal-metaethics-and-metasemantics-for-ai-alignment-2e72533cad6d
 5 
 6    Developed by June Ku in the set theoretic programming language setlX v2.2.0
 7    and interspersed with considerable non-technical philosophical commentary.
 8      https://randoom.org/Software/SetlX/ 
 9      https://download.randoom.org/setlX/pc/archive/setlX_v2-2-0.binary_only.zip
10        - contains among other things, tutorial.pdf introduction and reference
11 
12    I am very grateful to Howard Nye for our work together on metaethics as well 
13    as many helpful comments and conversations during much of this development.
14    I also appreciate helpful discussions with Mahendra Prasad, Jessica Taylor, 
15    Anna Salamon and Nate Soares.
16 */
17 load('lib.stlx');
18 load('causal_model.stlx');
19 load('causal_markov_model.stlx');
20 load('decision_algorithm.stlx');
21 
22 /* MetaEthical.AI Utility Function
23    ===============================
24    Takes a causal markov model of the world w and set of brains bs in w and 
25    returns a metaethical.ai utility function. 
26 
27    Abstract: We construct a fully technical ethical goal function for AI by 
28    directly tackling the philosophical problems of metaethics and mental 
29    content. To simplify our reduction of these philosophical challenges into 
30    "merely" engineering ones, we suppose that unlimited computation and a 
31    complete low-level causal model of the world and the adult human brains in 
32    it are available.
33 
34    Given such a model, the AI attributes beliefs and values to a brain in two 
35    stages. First, it identifies the syntax of a brainís mental content by 
36    selecting a decision algorithm which is i) isomorphic to the brainís causal 
37    processes and ii) best compresses its behavior while iii) maximizing charity.
38    The semantics of that content then consists first in sense data that 
39    primitively refer to their own occurrence and then in logical and causal 
40    structural combinations of such content.
41 
42    The resulting decision algorithm can capture how we decide what to *do*, but 
43    it can also identify the ethical factors that we seek to determine when we 
44    decide what to *value* or even how to *decide*. Unfolding the implications 
45    of those factors, we arrive at what we should do. All together, this allows 
46    us to imbue the AI with the necessary concepts to determine and do what we 
47    should want it to do.  
48 */
49 metaethical_ai_u := procedure(w, bs) {
50   /* Given a set of brain models, associate them with the decision algorithms
51      they implement.  */
52   d := { [b, decision_algorithm.implemented_by(b, {})] : b in bs };
53   /* Then map each brain to its rational self's values (understood extensionally
54      i.e. cashing out the meaning of their mental concepts in terms of the 
55      world events they refer to).  */
56   ext_rufs := { [b, d[b].ext_ruf(w,b)] : b in bs };
57 
58   state_space := w.cm().poss_states();
59   // Possible social choice utility functions
60   psc_us := choices({ {[state, r_i] : r_i in [1..#state_space] }
61                       : state in state_space });
62   sc_u := sort_set(psc_us, procedure(psc_u1, psc_u2) {
63                              return meai_psc_dist(psc_u1, ext_rufs) <
64                                     meai_psc_dist(psc_u2, ext_rufs);
65                            })[1];
66   return sc_u;
67 };
68 
69   meai_psc_dist := procedure(psc_u, ext_rufs) {
70     return +/ { ord_u_dist(card_u_to_ord_u(ext_ruf), card_u_to_ord_u(psc_u))**2
71                 : ext_ruf in ext_rufs };
72   };
73 
74 /* MIT License
75    -----------
76    Copyright 2019 June Ku 
77 
78    Permission is hereby granted, free of charge, to any person obtaining a copy 
79    of this software and associated documentation files (the "Software"), to deal
80    in the Software without restriction, including without limitation the rights 
81    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 
82    copies of the Software, and to permit persons to whom the Software is 
83    furnished to do so, subject to the following conditions:
84 
85    The above copyright notice and this permission notice shall be included in 
86    all copies or substantial portions of the Software.
87 
88    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 
89    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 
90    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 
91    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 
92    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 
93    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS 
94    IN THE SOFTWARE.
95 */
   1 class decision_algorithm(i, o, p, cp, u, m, m2, n, r, t_s, t_o,
   2                          cached_f, cached_lf, cached_ref) {
   3 // Definition and Structure of Decision Algorithms 
   4 // ===============================================
   5 
   6   i := i; /* Set of input variables (or technically names, ie strings)
   7                which are a subset of range(p)
   8              e.g. { 'P(I_1)', 'P(I_2)', 'P(I_3)' } 
   9              i.e. the agent's sense data. 
  10              Generally, p should include it, e.g. { ['I_1', 'P(I_1)'], ... }
  11   Technical footnote: There's some abuse of notation here. Each string like 
  12   'I_1' can depending on context be treated as:
  13     i) an individual constant term / object. A predicate expression in the 
  14        domain of p that has this string as an argument would be the 
  15        prototypical context illustrating this usage. Semantically, these would 
  16        refer to the brain events that implement them.  
  17     ii) a couple possibilities for this option: 
  18       a) a 0-ary predicate or sentential expression. It takes no inputs and 
  19          forms a proposition on its own. 
  20       b) shorthand for a 1-place predicate that takes the constant in ii and 
  21          simply asserts its existence, perhaps loosely inspired by Free Logic's
  22          E!. 
  23       In either case, its truth-maker/referent would be the same as option i's.
  24       Treating the string as a sentential expression or forming an expression 
  25       with it as the argument of a connective would indicate this usage.  
  26     All of these possibilities are expressions, or terms that could figure in an
  27     expression, i.e. the sort of thing that would be in the domain of p. They 
  28     should be distinguished from the set i, which would be a subset of the range
  29     of p. These would be names of a decision algorithm variable, e.g. parents() 
  30     of other variables. Assignments of value to this variable would be a 
  31     subjective probability 0-1. That would form a local/partial decision state,
  32     which could get mapped to by a local/partial implementation functor lf() 
  33     from a local brain state.  */
  34 
  35   o := o; /* Set of output variables
  36              e.g. { 'O_1', 'O_2', 'O_3' }
  37              i.e. the agent's motor output. 
  38              See also the technical note for i.  */
  39 
  40   /* Function from logical formulas (setlx functors/strings/terms) to 
  41        probability variables
  42      e.g. { ['A', 'P(A)'], ['B', 'P(B)'], [And('A', 'B'), 'P(A & B)'], 
  43             [BoxArrow('A', 'B'), 'P(A []-> B)'] }  */
  44   p := p;
  45 
  46   /* Function from ordered pairs [x, y] to conditional probability variables
  47      e.g. { [[x, y], 'P(x|y)'], ... } 
  48      or more concretely: { [['C', And('A', 'B')], 'P(C|A,B)'], ... }  */
  49   cp := cp;
  50 
  51   /* Function from logical formulas to utility variables 
  52        e.g. { ['A', 'U(A)'], ... } 
  53      u has the same structure as p above.
  54      See also accepted norms / higher-order utilities a few lines below. */
  55   u := u;
  56 
  57   /* Set of memory variables 
  58      This is a catch-all for other cognitive states that don't fit into the 
  59      other sets of variables yet may be important to model, e.g. intermediate 
  60      states in a calculation performing inference or decision. One can think of
  61      it as roughly analogous to a computer's RAM.  */
  62   m := m;
  63 
  64   /* Set of memory variables for t_o 
  65      Same idea as m above but these only mediate between p + cp + u and o. 
  66      This helps ensure it is a genuinely decision algorithmic explanation.  */
  67   m2 := m2;
  68 
  69   /* List of accepted norms / higher-order utilities
  70      [ n_1, n_2, ... ]
  71      each n_i == { ['e', e_n_i], ['m2', m2_n_i], ['o', o_n_i] }
  72        e_n_i is a function from logical formulas to evaluation variables
  73          e.g. { ['A', 'E_n_i(A)'], ... }
  74          They are basically utility function variables but since they are
  75          higher-order, i.e. they influence normative judgments or prescriptions
  76          rather than motor output, we call them evaluation variables instead
  77          when we want to emphasize that they don't directly govern action.
  78        m2_n_i is a set of memory variables that only moderate p + e_n_i to 
  79          o_n_i
  80        o_n_i is a set of output variables that causally influence various u or
  81          e_n_j. In simplified models, we can make o_n_i trivially cause u or 
  82          e_n_i states by just being identical to them 
  83 
  84      These are essentially decision criteria that govern how we change our
  85      preferences or other higher-order norms. The way in which they govern is 
  86      analogous and mathematically isomorphic to how first-order utility 
  87      functions or preferences govern action. There, beliefs that some motor
  88      output would cause some preferred state to be achieved will tend to cause
  89      such output and therefore the desired state. In the case of these accepted
  90      norms, beliefs that the output variables (which we can also understand as 
  91      normative judgments or prescriptions) will cause or constitute a change 
  92      in values (which is what the exprs of the domain of e_n_i would refer to) 
  93      will tend to cause such outputs and therefore changes in value. 
  94 
  95      It is very important to note that what makes these preferences / utilities 
  96      higher-order is their causal or functional roles and not the content of 
  97      the expressions associated with the evaluation variables. One can have, in 
  98      our sense, first-order preferences which prefer states involving other
  99      preferences. Those might be referentially "higher-order" in some sense but 
 100      they would still be first-order in our sense so long as they govern
 101      actions rather than normative judgments or prescriptions.
 102 
 103      See http://www.metaethical.ai/norm_descriptivism.pdf for more.
 104 
 105      One might object that introducing these higher-order utility functions
 106      invalidates what was supposed to be the utility function of the decision
 107      algorithm. Since a utility function is supposed to encompass any factors   
 108      explaining behavior, wouldn't it already include these allegedly
 109      higher-order factors? Strictly speaking, that may be correct, but that kind
 110      of utility function is not the most useful for our purposes.
 111      
 112      It may be better to think in terms of modeling a decision algorithm as a
 113      society of subagents. The first-order utility function then is the utility
 114      function of the subagent whose outputs are the decision algorithm's actions
 115      rather than prescriptions that directly influence other subagents' 
 116      utilities. When we speak of the utility function of the decision algorithm,
 117      this should just be understood as shorthand for this first-order utility 
 118      function.
 119   */
 120   n := n;
 121 
 122   // Internal states
 123   s := procedure() {
 124     if (n == {} || n == [] || n == om) {
 125       n_vs := {};
 126     } else {
 127       n_vs := +/ { range(n_i['e']) + n_i['m2'] + n_i['o'] : n_i in n };
 128     }
 129     return range(p) + range(cp) + range(u) + m + m2 + n_vs - i;
 130   };
 131 
 132   // Function from i + s() + o to ranges of values for those variables
 133   r := r;
 134 
 135   /* Higher order function from a function assigning values to all input and 
 136      internal states to a function assigning values to all next internal states 
 137        is -> s
 138   */
 139   t_s := t_s;
 140 
 141   /* Higher order function from a function assigning values to all probability 
 142      and utility variables as well as m2 to a function assigning values to 
 143      outputs 
 144        p u m2 -> o
 145   */
 146   t_o := t_o;
 147 
 148   /* State transitions
 149      Returns a higher-order function from functions assigning values to i + s 
 150      to functions assigning values to s + o. This is the intuitive combination
 151      of t_s and t_o. */
 152   t := procedure() {
 153     return { [    // a higher-order function which takes 
 154               is, // a function assigning values to i + s and returns a function
 155               t_s[is] + t_o[just_p_u_m2(is)]
 156                   // composed of the union of a func assigning values to s and 
 157                   // a func assigning values to o
 158              ] : is in domain(t_s) };
 159   };
 160 
 161   /* We explicitly cache f after computing it. This allows us to directly set
 162      its value when testing or when performing the brute force search would be
 163      prohibitively costly.  */
 164   cached_f := cached_f;
 165   cached_lf := cached_lf;
 166 
 167   /* ns as in the plural of n. 
 168      Takes 'e', 'm2' or 'o' and returns the corresponding variables from all n.
 169      e.g. ns('e') would return all e_n_i */
 170   ns := procedure(str) {
 171     if (str == 'e') {
 172       return +/ { range(n_i[str]) : n_i in n };
 173     } else {
 174       return +/ { n_i[str] : n_i in n };
 175     }
 176   };
 177 
 178   p_vars := procedure() {
 179     return range(p) + range(cp);
 180   };
 181 
 182   p_u_vars := procedure() {
 183     return p_vars() + range(u);
 184   };
 185 
 186   /* Takes a function (is) assigning values to input and internal state vars
 187      and returns its partial function with only vars in p, u and m2  */
 188   just_p_u_m2 := procedure(is) {
 189     return { [var, val] : [var, val] in is | var in p_u_vars() + m2 };
 190   };
 191 
 192 // Variable Relationships 
 193 // ----------------------
 194 
 195   /* Takes a variable and returns its parents, i.e. the variables that 
 196      are necessary to determine its next state in the state transition t(). */
 197   parents := procedure(y) {
 198     result := { z : z in (i + s()) | is_necessary_parent_of(z, y) };
 199     return result;
 200   };
 201 
 202     is_necessary_parent_of := procedure(z, y) {
 203       vs := {v : v in i + s() | v != z}; // all other variables but z
 204       f_vs := possible_states(vs, r);
 205       bool := exists(f_v in f_vs | vals_diff(f_v, y, z));
 206       return bool;
 207     };
 208 
 209     /* Check whether there exists an assignment of values to the other 
 210        variables such that there are still different values of y depending 
 211        on values of z.  */
 212     vals_diff := procedure(f_v, y, z) {
 213       vals := { t_s[ f_v + {[z,r_i]} ][y] : r_i in r[z] };
 214       return #vals > 1;
 215     };
 216 
 217   // Takes a variable and returns its children
 218   children := procedure(y) {
 219     return { v : v in s() + o | y in parents(v) };
 220   };
 221 
 222   /* Takes a variable v and a set of variables vs (with v usually in vs)
 223      and returns the first ancestors within any given branch of v's ancestors 
 224      that are not in vs */
 225   non_shared_ancestors := procedure(v, vs) {
 226     result := {};
 227     for (parent in parents(v)) {
 228       if (parent in vs) {
 229         result := result + non_shared_ancestors(parent, vs);
 230       } else {
 231         result := result + { parent };
 232       }
 233     }
 234     return result;
 235   };
 236 
 237   /* Takes a variable v and a set of variables vs (with v usually in vs)
 238      and returns the first descendants within any given branch of v's 
 239      descendants that are not in vs */
 240   non_shared_descendants := procedure(v, vs) {
 241     result := {};
 242     for (child in children(v)) {
 243       if (child in vs) {
 244         result := result + non_shared_descendants(child, vs);
 245       } else {
 246         result := result + { child };
 247       }
 248     }
 249     return result;
 250   };
 251 
 252 
 253 // Syntax and Implementation
 254 // =========================
 255 // See also better_explanation and the section starting with complexity.
 256 
 257   /* Function from a causal_markov_moodel of a brain to an implementation 
 258      mapping f which is a higher-order function from a function assigning 
 259      states to the brain's variables to a function assigning states to internal 
 260      state variables s() meeting certain criteria. 
 261   */
 262   f := procedure(b) {
 263     if (cached_f == om) {
 264       cached_f := compute_f(b);
 265     }
 266     return cached_f;
 267   };
 268 
 269   compute_f := cachedProcedure(b) {
 270     fs := { pf : pf in possible_fs(b) | commutes(pf, b) };
 271     /* Todo: Sort by complexity and return the simplest? 
 272        Or also factor in ambitiousness? Coverage of more states?  */
 273     if (#fs > 0) {
 274       return first(fs);
 275     } else {
 276       return {};
 277     }
 278   };
 279 
 280   /* Function from a causal_markov_model of a brain to the possible 
 281      implementation mappings f which map from a function assigning values to 
 282      the exogenous and endogenous variables of the causal_markov_model to 
 283      functions assigning values to the input, internal state, and output 
 284      variables. */
 285   possible_fs := procedure(b) {
 286     return function_space( b.poss_states(), poss_states() );
 287   };
 288 
 289   poss_states := cachedProcedure() {
 290     return possible_states(i + s() + o, r);
 291   };
 292 
 293   is_implemented_by := procedure(b) {
 294     return !(f(b) == {});
 295   };
 296 
 297   /* Implementation Functor
 298      ----------------------
 299      Check whether a possible implementation functor pf commutes. See f.
 300 
 301      This basically implements Chalmers' criterion for when some causal system 
 302      implements a computation: http://consc.net/papers/rock.html
 303 
 304      In category theory terms, we can say that f is an implementation functor.
 305      Given a category of brain states and the causal relations between them
 306      and another of decision states and cognitive state transitions between
 307      them, f would be a functor mapping i) brain states to decision states and
 308      ii) causal transitions to decision state transitions such that it doesn't
 309      matter whether one goes from a) a brain state to the brain state it causes 
 310      and then to the decision state f maps that brain state to or b) first from 
 311      the same original brain state to the decision state f says it implements 
 312      and then to the next decision state this decision algorithm transitions to.
 313      Either way, you would end up at the same final decision state. 
 314   */
 315   commutes := procedure(pf, b) {
 316     /*                 pf
 317       bs_{b.uv_t_i-1} ------> ds_iso_t_i-1
 318        |                       |
 319        | b.response            | restrict_domain
 320        |                       | t()
 321        v                       v
 322       bs_{b.v_ti} ----------> ds_so_t_i
 323            +u  pf  restrict_domain
 324 
 325       Roughly check that t( pf ) == pf( b.response ) 
 326     */
 327     return forall(bs in domain(pf) |
 328        t()[restrict_domain(pf[bs], i + s())] ==
 329        restrict_domain(pf[
 330          b.response(b.v, bs) + { [u_var, first(b.r[u_var])] : u_var in b.u }
 331        ], s() + o)
 332      );
 333   };
 334 
 335   /* Takes a decision state assigning values to i + s() + o and returns a
 336      function assigning values to just s() + o  */
 337   drop_i := procedure(ds) {
 338     return { [var, val] : [var, val] in ds | !(var in i) };
 339   };
 340 
 341   /* Local / partial f
 342      Returns a function mapping partial causal states to partial decision states
 343      which commutes, i.e. it maps partial causal states to the partial decision 
 344      state which is forced by the complete mapping f.  */
 345   lf := procedure(b) {
 346     if (cached_lf == om) {
 347       cached_lf := compute_lf(b);
 348     }
 349     return cached_lf;
 350   };
 351 
 352   compute_lf := cachedProcedure(b) {
 353     /* Possible local fs 
 354        Function space from union of partial causal states to the union of 
 355        partial decision states.  */
 356     plfs := function_space(
 357       +/ { partial_functions(cs) : cs in b.poss_states() }, // causal states
 358       +/ { partial_functions(ds) : ds in poss_states() } // decision states
 359     );
 360     lfs := [ plf : plf in plfs | lf_commutes(plf, b) ];
 361     lfs := sort_list(lfs, procedure(plf1, plf2) { return #plf1 > #plf2; });
 362     return lfs[1];
 363   };
 364 
 365   /* Function from possible local function plf and a causal model of brain b 
 366      to boolean value whether it commutes. */
 367   lf_commutes := procedure(plf, b) {
 368     return forall(pcs in domain(plf) |  // pcs = partial causal state
 369       /* The set of complete decision states that is mapped by f from each 
 370          complete causal state compatible with the partial causal state is 
 371          equal to the set of complete decision states compatible with 
 372          the partial decision states mapped by the possible lf from the partial
 373          causal state. */
 374       { f(b)[ccs] : ccs in compatible_complete_states(pcs, b.poss_states()) } ==
 375       compatible_complete_states(plf[pcs], poss_states())
 376     );
 377   };
 378 
 379   /* Local / partial f inverse
 380   /* Returns a "function" from local / partial decision states to a set of 
 381      local / partial causal states, namely the inverse of lf.  */
 382   lf_inv := cachedProcedure(b) {
 383     return { [ds, { cs : cs in domain(lf(b)) | lf(b)[cs] == ds} ] :
 384              ds in range(lf(b)) };
 385   };
 386 
 387   // Helper function. Same as lf_inv but time subscripts are added to the b vars
 388   lf_inv_at := procedure(b, t_i) {
 389     return { [ds, { assign_time(cs, t_i) : cs in css }]
 390              : [ds, css] in lf_inv(b) };
 391   };
 392 
 393   /* Takes an event in b.cm() and returns the local decision state that the 
 394      brain state implements.  */
 395   lf_cm := procedure(b, b_event) {
 396     return lf(b)[drop_time(b_event)];
 397   };
 398 
 399   // Returns the decision state variables implemented by a local brain state.
 400   lf_cm_vs := procedure(b, b_event) {
 401     d_event := lf_cm(b, b_event);
 402     if (d_event == om) {
 403       return {};
 404     } else {
 405       return domain(d_event);
 406     }
 407   };
 408 
 409   /* Similar to above but returns the formula in the domain(p) that gets 
 410      mapped to the p var.  */
 411   lf_cm_p_inv := procedure(b, b_event) {
 412     vs := lf_cm_vs(b, b_event);
 413     if (#vs == 1) {
 414       return p_inv(first(vs));
 415     } else {
 416       return om;
 417     }
 418   };
 419 
 420   plf_to_pf := procedure(b, plf) {
 421     /* For each decision state var, construct the function from complete 
 422        causal states to the state of that variable. */
 423     pf_var := {
 424       [var, +/ { { [ccs, pds]
 425                    : ccs in compatible_complete_states(pcs, b.poss_states())
 426                  } : [pcs, pds] in plf | first(first(pds)) == var }
 427       ] : var in i + s() + o
 428     };
 429     return { [bs, +/ { (pf_var[var])[bs] : var in i + s() + o }]
 430              : bs in b.poss_states() };
 431   };
 432 
 433 
 434   // This is mostly just checking for syntactic validity or well-formedness.
 435   // See also all_lte.
 436   is_valid := procedure() {
 437     return m2s_are_valid() && no_e_in_u();
 438   };
 439 
 440   /* Ensures that m2s merely moderate between p + u and o.
 441      I'm no longer sure if this is necessary since we measure and punish 
 442      attributing irrationality but it also doesn't seem too bad of a 
 443      simplifying assumption for now.  */
 444   m2s_are_valid := procedure() {
 445     return forall (m2_i in m2 |
 446       forall (v in non_shared_ancestors(m2_i, m2) | v in p_u_vars()) &&
 447       forall (v in non_shared_descendants(m2_i, m2) | v in o)
 448     ) && forall (n_i in n |
 449       forall (m2_j_n_i in n_i['m2'] |
 450         forall (v in non_shared_ancestors(m2_j_n_i, n_i['m2']) |
 451                 v in p_vars() + range(n_i['e']) ) &&
 452         forall (v in non_shared_descendants(m2_j_n_i, n_i['m2']) |
 453                 v in n_i['o'])
 454       )
 455     );
 456   };
 457 
 458   /* This helps ensure higher-order ambitiousness by preventing the smuggling
 459      of evaluation functions within u. For each p expr of the form x []-> y,
 460      make sure it does not take the form that x-ing would bring about y which 
 461      the agent values (but x is not in the normal outputs o) and this 
 462      recognition causes x-ing. If it is of that form, we want a decision 
 463      algorithm that includes that explicitly within n.  
 464      Todo: Generalize to no_e_in_e and cover more cases, e.g. indirect 
 465      reference or equivalent expressions or influencing as grandparents.  */
 466   no_e_in_u := procedure() {
 467     return forall(expr in domain(p) | fct(expr) != "BoxArrow" || !e_in_u(expr));
 468   };
 469 
 470     e_in_u := procedure(expr) {
 471       x := args(expr)[1]; // like an n[i]['o']
 472       y := args(expr)[2];
 473       return y in domain(u) && !(x in o) && p[expr] in parents(x);
 474     };
 475 
 476 // Semantics of Mental Content 
 477 // ===========================
 478 
 479   /* Takes a set of pairs [w, b] and returns a function from possible worlds 
 480      to extensions / referents. 
 481        See: https://plato.stanford.edu/entries/two-dimensional-semantics/ 
 482   */
 483   intension := procedure(poss_wbs) {
 484     return { [wb[1], ref(wb[1], wb[2])] : wb in poss_wbs };
 485   };
 486 
 487   /* Reference
 488      ---------
 489      Reference function from causal markov model of world w and brain b 
 490      to a function which takes a brain and returns its set of truth-makers if 
 491      any, i.e. a set of poss events in the world or equivalently, functions
 492      assigning values to vars in the world. 
 493 
 494      There are a few ways that reference is grounded.
 495 
 496        1. Input variable states refer to the brain states implementing them.
 497           You can think of their content as something like "this sense datum
 498           is occurring" or perhaps even as just the pure demonstrative "this".
 499           It seems intuitive that these states can represent themselves since
 500           everything trivially carries information about itself. 
 501 
 502        2. The meanings of logical connectives (e.g. And, Or, Implies) and the 
 503           causal connective BoxArrow are grounded by the axioms that define 
 504           them. These meanings are enforced by punishing any attributions of 
 505           these connectives into credences that violate their axioms. See 
 506           incoherence. Cf inferential / conceptual role semantics:
 507             https://plato.stanford.edu/entries/content-narrow/#rol  
 508 
 509        3. The rest can be understood as a kind of ramsification. Other posited 
 510           objects/events, predicates and relations refer to those things
 511           (if any) which play the roles posited for them. 
 512              https://en.wikipedia.org/wiki/Ramsey%E2%80%93Lewis_method
 513              http://www.jimpryor.net/teaching/courses/mind/notes/ramseylewis.html
 514           We soften the requirement of strict truth by scoring by squared error.
 515           This amounts to a principle of charity favoring interpretations
 516           that make more beliefs turn out to be true (while also adhering to
 517           the constraints above).  
 518   */
 519   ref := procedure(w, b) {
 520     if (cached_ref == om) {
 521       cached_ref := compute_ref(w, b);
 522     }
 523     return cached_ref;
 524   };
 525 
 526   compute_ref := cachedProcedure(w, b) {
 527     /* Base expressions are strings or predicates. We'll build composite 
 528        expressions together with corresponding references out of combinations of
 529        these and logical / causal connectives. */
 530     base_exprs := { [expr, e_arity] : [expr, e_arity] in str_exprs() + preds()
 531                                     | !(p[expr] in i) };
 532     pbers := poss_base_expr_refs(w, base_exprs);
 533     poss_atom_refs := { i_ref(w, b) + base_ref
 534                         : base_ref in poss_base_refs(w, b, pbers) };
 535     poss_refs := [ full_ref(w, b, poss_atom_ref)
 536                    : poss_atom_ref in poss_atom_refs ];
 537     return sort_list(poss_refs, less_sq_err)[1];
 538   };
 539 
 540     /* Input vars refer to their own brain-state implementors. */
 541     i_ref := procedure(w, b) {
 542       return {
 543         [b_event, compatible_complete_states(b_event, w.cm().poss_states())]
 544         : b_event in b.actual_sync_events()
 545         | #lf_cm(b, b_event) == 1 && lf_cm(b, b_event) in poss_i_events()
 546       };
 547     };
 548 
 549     str_exprs := procedure() {
 550       return { [expr, 0] : expr in domain(p)
 551                          | !(p[expr] in i) && isString(expr) };
 552     };
 553 
 554     preds := procedure() {
 555       return { [fct(expr), arity(expr)]
 556                : expr in domain(p)
 557                | isTerm(expr) && !(fct(expr) in connectives()) };
 558     };
 559 
 560     // Set of choices of poss assignments of world w states/sets/tuples to expr
 561     poss_base_expr_refs := procedure(w, base_exprs) {
 562       return choices({ {expr} >< range_for_arity(e_arity, w.cm().poss_states())
 563                        : [expr, e_arity] in base_exprs });
 564     };
 565 
 566     /* Use poss_base_expr_refs to construct poss reference functions assigning
 567        reference to brain events rather than exprs.  */
 568     poss_base_refs := procedure(w, b, pbers) {
 569       return {
 570         +/ {
 571           // For each brain event that represents the given expr
 572           // Assign to it the reference of that expr according to pber
 573           { [b_event, pber[expr]] : b_event in b.actual_sync_events()
 574                                   | lf_cm_p_inv(b, b_event) == expr }
 575           // Do that for each base expr and collect them into a single function
 576           : expr in domain(pber) }
 577         // for each poss base expr reference func
 578         : pber in pbers
 579       };
 580     };
 581 
 582     // Full Reference
 583     // Much of the below builds on fairly standard first-order structures or 
 584     // models, e.g., https://en.wikipedia.org/wiki/First-order_logic#Semantics 
 585     full_ref := procedure(w, b, atom_ref) {
 586       /* Full ref maps each brain event which implements a credence to the 
 587          result of reducing the reference of the corresponding expression. */
 588       return { [b_event, rr(
 589                             lf_cm_p_inv(b, b_event), // expr it is credence in 
 590                             [ w, b, atom_ref,
 591                               // vars are all the same time so just take 1st
 592                               var_time(first(first(b_event))),
 593                               {} ] // empty context
 594                          ) ]
 595                : b_event in actual_p_events(b) };
 596     }; // end full_ref
 597 
 598       /* Reduce reference, takes an expression and an array of further
 599          parameters consisting of a world, brain, atomic reference function, 
 600          time and context and returns a set of truth-making w events or the set 
 601          of w events the expression is true in. The set of truthmakers is our 
 602          way of enforcing a canonical boolean algebraic form. The different 
 603          elements act as disjuncts of events and combining assignments of vals 
 604          to vars into a map acts as conjunction. It may seem counterintuitive 
 605          and take some getting used to but we use this canonical form even when
 606          there there is only one disjunct. */
 607       rr := procedure(expr, arr) {
 608         w := arr[1]; b := arr[2]; atom_ref := arr[3]; ti := arr[4];
 609         context := arr[5];  // context is for binding of variables, see Exists
 610         if (isTerm(expr)) {
 611           if (#args(expr) >= 1) { a1 := args(expr)[1]; }
 612           if (#args(expr) >= 2) { a2 := args(expr)[2]; }
 613         }
 614         switch {
 615           case expr in domain(context) :
 616             return context[expr];
 617           case p[expr] in i :
 618             return compatible_complete_states(
 619                      first({ b_event : b_event in pow(b.a[ti])
 620                                      | lf_cm_p_inv(b, b_event) == expr }),
 621                      w.cm().poss_states()
 622                    );
 623           case isString(expr) :
 624             return atom_ref[expr];
 625           case isSet(expr) :
 626             /* This shouldn't happen but just in case, return the set.
 627                Maybe throw an error instead? Or it's probably just rr(rr(.. */
 628             print('Warning: Expected an expr but received a set: ' + expr);
 629             return expr;
 630           case fct(expr) == 'At_t' :
 631             if (ti + a2 >= 0) { // a2 should already be negative
 632               return rr(a1, [w, b, atom_ref, ti + a2, context]);
 633             } else {
 634               return {};
 635             }
 636           case fct(expr) == 'And' :
 637             return rr(a1, arr) * rr(a2, arr);
 638           case fct(expr) == 'Or' :
 639             return rr(a1, arr) + rr(a2, arr);
 640           case fct(expr) == 'Not' :
 641             if (isSet(a1)) { // Set of states
 642               return w.cm().poss_states() - a1;
 643             } else {
 644               // String or functor, i.e. predicate or relation application
 645               // [~A] = [~[A]]
 646               return rr(Not(rr(a1, arr)), arr);
 647             }
 648           case fct(expr) == 'Implies' :
 649             return rr( Or(Not(a1), a2), arr );
 650           case fct(expr) == 'Forall' :
 651             // Vx Px = ~Ex ~Px
 652             return rr(Not(Exists(a1, Not(a2))), arr);
 653           case fct(expr) == 'Exists' :
 654             // Set of truth-makers for a2 when any {w_event} is subbed for a1
 655             return +/ { rr(a2, [ w, b, atom_ref, ti,
 656                                  context + { [a1, {w_event}] } ])
 657                         : w_event in range_for_arity(0, w.cm().poss_states()) };
 658           case fct(expr) == 'BoxArrow' : // A []-> B, subjunctive conditional
 659             // Antecedent
 660             if (isSet(a1)) {
 661               ant := a1;
 662             } else {
 663               ant := rr(a1, arr);
 664             }
 665             // Consequent
 666             if (isSet(a2)) {
 667               cons := a2;
 668             } else {
 669               cons := rr(a2, arr);
 670             }
 671             if (#ant > 1) { // Antecedent has disjuncts.
 672               // This or that causes the effect if this does or that does.
 673               return +/ { rr(BoxArrow(disjunct, cons), arr) : disjunct in ant };
 674             } else if (#ant == 1) {
 675               if (#cons > 1) {
 676                 // A thing causes this or that if it causes this or causes that.
 677                 return +/ { rr(BoxArrow(ant, disjunct), arr) : disjunct in cons };
 678               } else if (#cons == 1) {
 679                 if ( w.cm().response( domain(first(cons)), first(ant) ) ==
 680                      first(cons) ) {
 681                   /* If this aspect of the causal model is true, it's true in
 682                      all states. Fire causes smoke can be true even in a     
 683                      possible world where all fire is prevented. A more 
 684                      expressive language could substitute a more meaningful 
 685                      truthmaker here. */
 686                   return w.cm().poss_states();
 687                 } else {
 688                   return {};
 689                 }
 690               } else {
 691                 // Trivially true because everything causes the empty effect? 
 692                 return w.cm().poss_states();
 693               }
 694             } else {
 695               return {}; // Nothing is caused by nothingness?
 696             }
 697             return om;
 698           default : // Predicate or Relation
 699             tuple := [rr(arg, arr) : arg in args(expr)];
 700             if (tuple in atom_ref[fct(expr)]) {
 701               // If a predication is true, it's true in all states.
 702               return w.cm().poss_states();
 703             } else {
 704               return {};
 705             }
 706         } // end switch
 707       }; // end rr
 708 
 709     less_sq_err := procedure(w, b, ref_1, ref_2) {
 710       return sq_err(w, b, ref_1) < sq_err(w, b, ref_2);
 711     };
 712 
 713     sq_err := procedure(w, b, ref_i) {
 714       return +/ {
 715         ( true_p(w, ref_i, b_event) - first(range(lf_cm(b, b_event))) ) ** 2
 716         : b_event in actual_p_events(b)
 717       };
 718     };
 719 
 720     true_p := procedure(w, ref_i, b_event) {
 721       if (w.actual_state() in ref_i[b_event]) {
 722         return 1;
 723       } else {
 724         return 0;
 725       }
 726     };
 727 
 728   // Takes an expression (in a world, brain and time) and returns its referent.
 729   ref_expr := procedure(w, b, t_i, expr) {
 730     p_var := p[expr];
 731     ds := f(b)[drop_time(b.a[t_i])]; // decision state
 732     p_event := { [p_var, ds[p_var]] };
 733     b_event := first( lf_inv_at(b, t_i)[p_event] );
 734     return ref(w, b)[b_event];
 735   };
 736 
 737   // Returns the brain events that implement a placement of probabibility.
 738   actual_p_events := procedure(b) {
 739     return { b_event : b_event in b.actual_sync_events()
 740              |      #lf_cm_vs(b, b_event) == 1 &&         // just one var
 741                first(lf_cm_vs(b, b_event)) in range(p) }; // and it's a p_var
 742   };
 743 
 744   /* Function from causal_markov_model b to a function from times to 
 745      probability states, i.e. assignments of values to range(p) + range(cp).
 746      At each time, use f(b) to see what decision state it maps the 
 747      brain state to and collect the states for range(p) + range(cp). 
 748 
 749      Actual world b.a is a function from times to assignments of values to 
 750      variables at that time.
 751         e.g. b.a := { [1, assign_time({ [y, true] : y in b.uv() }, 1)], ...} */
 752   prob_states := procedure(b) {
 753     return { [ti, prob_states_t(b, ti) ] : ti in [1..b.n] };
 754   };
 755 
 756   /* Probability States at Time ti 
 757      Takes a brain b at a time ti and returns the decision state that is 
 758      implemented by that brain state but restricted to just the probability 
 759      variables.  */
 760   prob_states_t := procedure(b, ti) {
 761     bs := { [var, val] : [var, val] in drop_time(b.a[ti]) | var in b.uv() };
 762     ds := f(b)[bs];
 763     return { [var, val] : [var, val] in ds | var in p_vars() };
 764   };
 765 
 766 
 767 // Utility / Evaluation Functions
 768 // ==============================
 769 
 770   // Collect all utility/evaluation function vars with their output vars
 771   eo := procedure() {
 772     return {[u, o]} + { [n_i['e'], n_i['o']] : n_i in n };
 773   };
 774 
 775   /* Instrumental irrationality, including of evaluation as well as utility 
 776      functions. */
 777   instr_irrat := procedure(w, b) {
 778     return +/ { e_dist(w, b, ti) : ti in [1..(b.n - 1)] } / (b.n - 1);
 779   };
 780 
 781     /* Takes a brain b and time ti and returns the distance between the actual 
 782        outputs and the optimal outputs according to the utility / evaluation 
 783        functions and []-> beliefs.  */
 784     e_dist := procedure(w, b, ti) {
 785       /* These are all weighted equally right now, but we may want to 
 786          consider different weights. Should first-order utility be treated 
 787          differently? What about highest order? Does complexity matter?  */
 788       return +/ { e_i_dist(w, b, ti, e_i, o_i) : [e_i, o_i] in eo() } / #eo();
 789     }; // end e_dist
 790 
 791     // Similar to above but for a specific util/eval function and output vars
 792     e_i_dist := procedure(w, b, ti, e_i, o_i) {
 793       poss_o_i_states := possible_states(o_i, r);
 794       actual_d_st0 := f(b)[ drop_time(b.a[ti]  ) ];
 795       actual_d_st1 := f(b)[ drop_time(b.a[ti+1]) ];
 796       actual_o_i_state := { [o_i_j, actual_d_st1[o_i_j]] : o_i_j in o_i };
 797       sorted_o_i_states := sort_list( [x : x in poss_o_i_states],
 798         procedure(x, y) {
 799           return exp_u(w, b, ti, actual_d_st0, e_i, x) >
 800                  exp_u(w, b, ti, actual_d_st0, e_i, y);
 801         } );
 802       ix := index_of(actual_o_i_state, sorted_o_i_states);
 803       /* Cardinal measures seem too easy to abuse so we convert them to an 
 804          ordinal measure. */
 805       return ix / #sorted_o_i_states;
 806     }; // end e_i_dist
 807 
 808     // Expected Utility
 809     exp_u := procedure(w, b, ti, actual_d_st, e_i, o_i_state) {
 810       /* Using our reference function, find an expression that refers to 
 811          the given output state. */
 812       o_expr_b_events := {
 813         b_event : b_event in actual_p_events(b)
 814                 | var_time(first(b_event)) == ti &&
 815                   ref(w, b)[b_event] != om &&
 816                   exists(w_state in ref(w, b)[b_event] |
 817                     exists(b_ev in lf_inv(b)[o_i_state] |
 818                       restricted_eq(w_state, assign_time(b_ev, ti+1))
 819                     )
 820                   )
 821       };
 822       assert(o_expr_b_events != {}, "No o_expr_b_events found");
 823       o_expr := lf_cm_p_inv(b, first(o_expr_b_events));
 824       // Todo: Handle multiple exprs referring to o_i_state
 825       /* Is the above too indirect and complicated? There may be worries of 
 826          not capturing the correct mode of presentation? Could we "hardcode"
 827          something into the language instead? Like with input vars or
 828          memories using At_t. Or the simplest may be to allow it to overlap
 829          with input vars. That is actually still allowed by the above. 
 830            Is this not so easy though since it needs to refer to a *future* o 
 831          state?
 832 
 833          Add up the expectations of utility. Here we are interpreting our   
 834          utility vars as additive values. U(x & y) = U(x) + U(y). If, 
 835          however, there are exceptions, they can still be encoded by 
 836          having a util var associated with the conjunction. Then, if U(x)
 837          and U(y) are still defined, you may need to be careful to encode 
 838          only the intended difference from additive in the conjunction 
 839          rather than doubly counting its full value. 
 840            Here, we are assuming the P(o []-> u) beliefs reflect the agent's
 841          judgment of overall utility. If the agent is poor at aggregating 
 842          utilities, their final judgment might come apart from these 
 843          credences. There may also be issues with the timing of these
 844          and any aggregated judgments. While this may require much 
 845          technical work to thoroughly address, we hope it is at least 
 846          conceptually and philosophical clear what is to be done. */
 847       return +/ {  actual_d_st[ p[BoxArrow(o_expr, u_expr)] ] // P(o []-> u)
 848                  * actual_d_st[e_i_j] : [u_expr, e_i_j] in e_i }; // * u
 849     }; // end exp_u
 850 
 851 // Probabilistic Coherence 
 852 // =======================
 853 
 854   incoherence := procedure(b) {
 855     /* Synchronic irrationality
 856        Avg of distances to the closest coherent probability distribution 
 857        at each time normalized by max distance. */
 858     sync_irrat := avg({ p_dist(b, ti) / max_p_dist(b, ti) : ti in [1..b.n] });
 859     /* Diachronic irrationality
 860        Probabilities either stay the same or get updated to be close to what
 861        a (locally) ideally rational inference would result in.  */
 862     dia_irrat := 1 - p_conn_continuity([ f(b)[drop_time(b.a[ti])]
 863                                          : ti in [1..b.n] ]);
 864     return avg([sync_irrat, dia_irrat]);
 865   };
 866 
 867   /* Maximum Probability Distance 
 868      Find the probability distribution that is the farthest distance from the 
 869      given brain's probability distribution and return that distance.
 870      This is used to normalize probability distance scores.  */
 871   max_p_dist := procedure(b, ti) {
 872     poss_p_states := [poss_p : poss_p in possible_states(p_vars(), r)];
 873     p_states := prob_states_t(b, ti);
 874     farthest_pp := sort_list(poss_p_states, procedure(pp1,pp2) {
 875                                return prob_distance(p, cp, p_states, pp1) >
 876                                       prob_distance(p, cp, p_states, pp2);
 877                              })[1];
 878     return prob_distance(p, cp, p_states, farthest_pp);
 879   };
 880 
 881   /* Returns the set of possible causal models formed from the base variables
 882      in the domain of p. */
 883   poss_cms := procedure() {
 884     base_vars := [ p[expr] : expr in domain(p) |
 885                    isString(expr) || (isTerm(expr) && fct(expr) == 'At_t') ];
 886     vs := { v : v in base_vars };
 887     orderings := permutations(base_vars);
 888     return +/ { poss_cms_from_l(l, vs) : l in orderings };
 889   };
 890 
 891     poss_cms_from_l := procedure(l, vs) { // l is an ordering (list)
 892       // Each var's poss parents are sets containing the preceding vars in l.
 893       poss_v_parents := { [l[j], pow({ v : v in l[1..j-1] })] : j in [2..#l] };
 894       poss_v_parents[l[1]] := {{}};
 895       poss_parents := possible_states(vs, poss_v_parents);
 896 
 897       return +/ { poss_cms_from_p_pars(vs, poss_parents, p_pars)
 898                   : p_pars in poss_parents };
 899     };
 900 
 901     /* Given a set of parenthood relations, construct its set of possible 
 902        causal models. For each variable, construct the set of possible 
 903        functions from its parents. */
 904     poss_cms_from_p_pars := procedure(vs, poss_parents, p_pars) {
 905       cm_u := { v : v in vs | p_pars[v] == {} };
 906       cm_v := vs - cm_u;
 907       cm_r := { [v, {true, false}] : v in vs };
 908       poss_v_cm_fs := {
 909         [v, function_space(possible_states(p_pars[v], cm_r), {true, false})]
 910         : v in cm_v };
 911       poss_cm_fs := possible_states(cm_v, poss_v_cm_fs);
 912       return { causal_model(cm_u, cm_v, cm_r, cm_f, om) : cm_f in poss_cm_fs };
 913     };
 914 
 915   /* Takes a causal_model and returns a set of coherent probability 
 916      distributions compatible with the cm. Each distribution is given by
 917      a pair of assignments of values to p and cp.  */
 918   coh_pcms := procedure(cm) {
 919     // Coherent probability distributions for exogenous variables u
 920     coh_u_probs := { p_u : p_u in poss_p_us(cm.u)
 921                          | p_u_is_coherent(p_u, cm.u) };
 922     return { [probs(cm, p_u) + subj_conds(cm, p_u), cond_probs(cm, p_u)]
 923              : p_u in coh_u_probs };
 924   };
 925 
 926     poss_p_us := procedure(us) {
 927       return choices({ {state} >< p_range() : state in state_space(us) });
 928     };
 929 
 930       p_range := procedure() {
 931         return r[first(range(p))]; // assume all probs share the same range
 932       };
 933 
 934       /* For each prob variable, choose either it being true or false */
 935       state_space := procedure(vs) {
 936         return choices({ {[v, true], [v, false]} : v in vs });
 937       };
 938 
 939     p_u_is_coherent := procedure(p_u, base_vars) {
 940       poss_exprs := da().add_up_to(base_vars, 2 ** #base_vars);
 941       poss_probs := possible_states(poss_exprs,
 942                                     { [expr, p_range()] : expr in poss_exprs });
 943       poss_ord_pairs := poss_exprs >< poss_exprs;
 944       poss_cond_probs := possible_states(poss_ord_pairs,
 945                                 { [pair, p_range()] : pair in poss_ord_pairs });
 946       poss_p_cps := poss_probs >< poss_cond_probs;
 947       return exists(p_cp in poss_p_cps
 948                     | is_superset_of(p_cp, p_u) && superset_is_coherent(p_cp) );
 949     };
 950 
 951       is_superset_of := procedure(p_cp, p_u) {
 952         return { [p_var, p_val] : [p_var, p_val] in p_u | p_val != om }
 953                < first(p_cp);
 954       };
 955 
 956       superset_is_coherent := procedure(p_cp) {
 957         d := decision_algorithm.new();
 958         d.p := { [expr, 'P(' + str(expr) + ')']
 959                  : [expr, p_val] in first(p_cp) };
 960         d.cp := { [l, 'P(' + str(l[1]) + '|' + str(l[2]) + ')']
 961                   : [l, p_val] in last(p_cp) };
 962         pp := { ['P(' + str(expr) + ')',                   p_val]
 963                 : [expr, p_val] in first(p_cp) } +
 964               { ['P(' + str(l[1]) + '|' + str(l[2]) + ')', p_val]
 965                 : [l, p_val] in last(p_cp) };
 966         return d.is_coherent(pp);
 967       };
 968 
 969     /* Given p_u, a prob distr over states of the exogenous variables u, use 
 970        the causal_model to calculate the probabilities of endogenous vars vs. */
 971     probs := procedure(cm, p_u) {
 972       // For each state of exogenous vars u, use the causal model to fill in 
 973       // what the state of endogenous vars v would be. 
 974       p_uv := { [u_state + cm.response(cm.v, u_state), p_val]
 975                 : [u_state, p_val] in p_u };
 976       // For each poss event of the causal model, turn the event into an expr
 977       // and assign it probability equal to the sum of the probabilities of 
 978       // uv states in which the event is realized.
 979       return { [ func_to_expr(g_vs),
 980                  +/ { p_val : [uv_state, p_val] in p_uv | g_vs < uv_state }]
 981                : g_vs in cm.poss_events() };
 982     };
 983 
 984     /* Next we calculate the conditional probabilities for a given p_u */
 985     cond_probs := procedure(cm, p_u) {
 986       return +/ { calc_cond_prob(g, h) : [g,h] in cm.poss_events() ** 2 };
 987     };
 988 
 989       calc_cond_prob := procedure(cm, p_u, g, h) {
 990         if (compatible_values(g, h)) {
 991           a := func_to_expr(g);
 992           b := func_to_expr(h);
 993           c := func_to_expr(g + h);
 994           /* Using P(A|B) == P(And(A,B)) / P(B)  (Kolmogorov definition) */
 995           return { [ cp[[a,b]], probs(cm, p_u)[c] / probs(cm, p_u)[b] ] };
 996         } else {
 997           return 0;
 998         }
 999       };
1000 
1001     // Just set to 100% if the subjunctive conditional is true in cm
1002     subj_conds := procedure(cm, p_u) {
1003       return { [p[BoxArrow(func_to_expr(g), func_to_expr(h))], 1.0]
1004                 : [g,h] in cm.poss_events() ** 2
1005                 | cm.response(domain(h), g) == h };
1006     };
1007 
1008   // Note that the function it returns assigns true/false to prob vars and not
1009   // atomic expressions. You may want expr_to_events instead.
1010   expr_to_p_funcs := procedure(expr, dom) {
1011     return { { [p[term], bool] : [term, bool] in event }
1012              : event in expr_to_events(expr, dom) };
1013   };
1014 
1015   /* Coherent Causal Model Probabilities
1016      Coherent probability distributions over p and cp, including the 
1017      subjunctive conditionals, which can be seen as pieces of causal models. 
1018      We start with probability distributions over all causal models (from
1019      the set of base variables in p). Collecting the coherent distributions 
1020      over p and cp within each model, we then distribute the probability of
1021      the causal model to those distributions and sum them all up again. */
1022   coh_cm_ps := procedure() {
1023     // Possible probabilities of causal models
1024     poss_p_cms := choices({ {poss_cm} >< p_range : poss_cm in poss_cms() });
1025     /* The probabilities within any distribution must sum to 1. 
1026        (Each causal model considered is global and includes all base vars so 
1027        there is no worries of overlapping by one cm being a "superset" of 
1028        another.) */
1029     coh_p_cms := { p_cms : p_cms in poss_p_cms | +/ range(p_cms) == 1.0 };
1030 
1031     // Collect from each coherent probability distribution over causal models
1032     // the coherent probability distributions over p and cp that are allowed.
1033     return +/ { sum_p_cms(p_cms) : p_cms in coh_p_cms };
1034 
1035     /* Todo: Double-check whether it's necessary to expand the probability 
1036        distributions to non-canonical expressions in p? Now that prob_distance
1037        handles logical equivalences, it seems unnecessary.  */
1038   };
1039 
1040     /* Given a coherent probability distribution over causal models... */
1041     sum_p_cms := procedure(p_cms) {
1042       /* For each causal model and its probability, collect the possible
1043          choices of a coherent probability distribution from each, where
1044          those probability distributions have had the probability of the 
1045          causal model distributed over them.  */
1046       sets_of_p_cps := choices({
1047         { distr_p_cm(p_cp, p_cm) : p_cp in coh_pcms(cm) }
1048         : [cm, p_cm] in p_cms
1049       });
1050 
1051       /* ... For each choice of p_cps from each causal model, sum them up. */
1052       return { foldr(add_p_distr, {}, [p_cp : p_cp in p_cps])
1053                : p_cps in sets_of_p_cps };
1054     };
1055 
1056       /* Distribute probability of causal model.
1057          Takes p_cp, a set of pairs of coherent probability distributions over p
1058          and cp, and a probability of the causal model p_cm and then 
1059          distributes the probability of the causal model so that the probability
1060          of each expression is reduced by multiplying the p_cm. */
1061       distr_p_cm := procedure(p_cp, p_cm) {
1062         // p_cp == [cm_p, cm_cp]
1063         cm_p := p_cp[1];
1064         cm_cp := p_cp[2];
1065         return [{ [p_var, p_val * p_cm] : [p_var, p_val] in cm_p },
1066                 { [cp_l, cp_val * p_cm] : [cp_l, cp_val] in cm_cp }];
1067       };
1068 
1069       /* Given two pairs of probability distributions over p_vars
1070          add up the probabilities of expressions. */
1071       add_p_distr := procedure(p_cp1, p_cp2) {
1072         p1  := p_cp1[1];
1073         cp1 := p_cp1[2];
1074         p2  := p_cp2[1];
1075         cp2 := p_cp2[2];
1076         return [sum_p(p1, p2), sum_p(cp1, cp2)];
1077       };
1078 
1079         sum_p := procedure(p1, p2) {
1080           /* If one has a variable where the other does not have any variable
1081              with an equivalent expression, just go with the probability 
1082              defined by the one. */
1083           return { [v1, p1[v1]] : v1 in domain(p1)
1084                    | !exists(v2 in domain(p2) | equiv_p_expr(p, cp, v1, v2)) } +
1085                  { [v2, p2[v2]] : v2 in domain(p2)
1086                    | !exists(v1 in domain(p1) | equiv_p_expr(p, cp, v1, v2)) } +
1087           /* Otherwise, for any variables that have equivalent expressions,
1088              add up the probabilities for those variables. (Since they've 
1089              already been reduced in distr_p_cm, we don't need to average 
1090              the sum.)  */
1091                  +/ { { [v1, p1[v1] + p2[v2]],[v2, p1[v1] + p2[v2]]  }
1092                       : [v1, v2] in domain(p1) >< domain(p2)
1093                       | equiv_p_expr(p, cp, v1, v2) };
1094         };
1095 
1096           equiv_p_expr := procedure(p, cp, v1, v2) {
1097             e1 := p_cp_inv(p, cp, v1);
1098             e2 := p_cp_inv(p, cp, v2);
1099             return equiv_expr(e1, e2, all_atomic_formulas(domain(p)));
1100           };
1101 
1102   /* Takes a brain and time and finds the closest perfectly coherent 
1103      probability distribution and returns the distance to it. */
1104   p_dist := procedure(b, ti) {
1105     p_states := prob_states_t(b, ti);
1106     return prob_distance(p, cp, sorted_rat(p_states)[1], p_states);
1107   };
1108 
1109   /* Takes a prob state and returns a list of coherent probability 
1110      distributions sorted by how close they are p_states, closer ones first. */
1111   sorted_rat := procedure(p_states) {
1112     return sort_list( [ coh_p + coh_cp : [coh_p, coh_cp] in coh_cm_ps() ],
1113       /* A procedure measuring which of two probability states is closer to 
1114          the original one given.  */
1115       procedure(pp1, pp2) {
1116         return prob_distance(p, cp, pp1, p_states) <
1117                prob_distance(p, cp, pp2, p_states);
1118       }
1119     );
1120   };
1121 
1122   /* Takes a p or cp var in range(p) + range(cp) and returns its preimage 
1123      formula(s) */
1124   p_inv := procedure(v) {
1125     if (v in range(p)) {
1126       return inv(p)[v];
1127     } else {
1128       return inv(cp)[v];
1129     }
1130   };
1131 
1132   // Same as p_inv but with p and cp explicitly passed in due to setlx bug
1133   p_cp_inv := procedure(p, cp, v) {
1134     if (v in range(p)) {
1135       return inv(p)[v];
1136     } else {
1137       return inv(cp)[v];
1138     }
1139   };
1140 
1141   // Same as p_inv but if in cp, it returns the first expr instead of the pair
1142   p_var_inv := procedure(v) {
1143     if (v in range(p)) {
1144       return inv(p)[v];
1145     } else {
1146       return first(inv(cp)[v]);
1147     }
1148   };
1149 
1150 
1151   /* Obeys the axioms of probability.
1152        See especially p 45 (or 64 in pdf) of E.T. Jaynes, Probability Theory: 
1153        The Logic of Science. https://bayes.wustl.edu/etj/prob/book.pdf  */
1154   old_is_coherent := procedure(pp) {
1155     /* For easier handling, turn domain(p) into a datatype like domain(cp) 
1156        so that e.g. P(A) looks more like P(A|om) */
1157     force_list := procedure(v) {
1158       if (isList(p_inv(v))) {
1159         return p_inv(v);
1160       } else {
1161         return [p_inv(v), om];
1162       }
1163     };
1164     pp := { [force_list(v), num] : [v,num] in pp };
1165     obeys_axiom_3 := procedure(l) {
1166       if (isTerm(l[1]) && fct(l[1]) == "And") { // [And(a,b), c]
1167         [a,b] := args(l[1]);
1168         c := l[2];
1169         // P(a,b|c) = P(a|b,c) * P(b|c)
1170         return pp[l] == pp[ [a, And(b,c)] ] * pp[ [b,c] ];
1171       } else {
1172         return true;
1173       }
1174     };
1175     return forall(l in domain(pp) | pp[l] >= 0) &&
1176            // P(A|x) + P(~A|x) = 1
1177            // Todo: Simplify double negation to avoid infinite NOT() regression?
1178            forall(l in domain(pp) | pp[l] == 1 - pp[ [Not(l[1]), l[2]] ]) &&
1179            forall(l in domain(pp) | obeys_axiom_3(l));
1180   }; // end old_is_coherent
1181 
1182   /* From Lasky, Kathryn. "Axiomatic First-Order Probability" 
1183      http://ceur-ws.org/Vol-527/paper5.pdf 
1184 
1185       Laskey axioms
1186         P1. probabilities are between 0 and 1
1187         P2. probabilities in tautologies are 1
1188               or just logical truths
1189         P3. if probability in two statements both being true are 0, then 
1190               probability of one or the other is equal to adding probability 
1191               in each
1192         P4. Bayesian conditioning
1193         P5. interchangeability of logically equivalent statements
1194 
1195      Todo: Double-check if it's necessary to handle undefined/om probability 
1196            values. It seems p_u_is_coherent checks its coherence by seeing if
1197            a superset is perfectly coherent so I think that strategy should 
1198            work. But there may be some slight technical issues making sure 
1199            we're allowing or disallowing om in the right places there.
1200   */
1201   is_coherent := procedure(pp) {
1202     pp := { [force_list(v), num] : [v,num] in pp };
1203     return forall(l in domain(pp) |
1204       // P1: Probabilities are between 0 and 1.
1205       0 <= pp[l] && pp[l] <= 1 &&
1206       // P2: Logical truths have probability 1.
1207       (!is_logical_truth(l[1]) || pp[l] == 1) &&
1208       forall(l2 in domain(pp) |
1209         ( l[2] != l2[2] || (
1210           // z is shared in P(_|z) && P(_|z)
1211 
1212           // P3: If P(x^y|z) = 0, then P(x v y|z) = P(x|z) + P(y|z)
1213           ( pp[[And(l[1], l2[1]), l[2]]] != 0 ||
1214             pp[[ Or(l[1], l2[1]), l[2]]] == pp[l] + pp[l2] ) &&
1215 
1216           // P4: Bayesian conditioning. P(x^y|z) = P(y|z) * P(x|y,z)  
1217           ( fct(l[1]) != "And" || !(l2[1] in args(l[1])) ||
1218             forall(l3 in domain(p) |
1219               !(l3[1] in args(l[1])) || l3[1] == l2[1] ||
1220               l3[2] != And(l2[1], l[2]) || pp[l] == pp[l2] * pp[l3]
1221             ) ) &&
1222 
1223           // P5: Interchangeability of equivalent expressions. Part 1
1224           //     If x<->y, then Vz P(x|z) = P(y|z) 
1225           ( !equiv_expr(l[1], l2[1], atomic_formulas(domain(p))) ||
1226             pp[l] == pp[l2] )
1227         ) ) /* end || */ && (
1228           // P5 Part 2: If x <-> y, then Vz P(z|x) = P(z|y) 
1229           !equiv_expr(l[2], l2[2], atomic_formulas(domain(p))) ||
1230           l[1] != l2[1] || pp[l] == pp[l2]
1231         )
1232       ) // end forall l2
1233     ); // end forall l
1234   }; // end is_coherent
1235 
1236     /* For easier handling, turn domain(p) into a datatype like domain(cp) 
1237        so that e.g. P(A) looks more like P(A|om) */
1238     force_list := procedure(v) {
1239       if (isList(p_inv(v))) {
1240         return p_inv(v);
1241       } else {
1242         return [p_inv(v), om];
1243       }
1244     };
1245 
1246   /* Takes two functions mapping p & cp vars to probabilities and computes
1247      the difference. See Staffel, Julia. Measuring the Overall Incoherence of 
1248      Credence Functions. Synthese, 2015. DOI 10.1007/s11229-014-0640-x 
1249 
1250      Technical Note: Takes p and cp (which should just be the usual this.p and 
1251      this.cp) but needs to be passed explicitly because of a setlx bug: The 
1252      second and subsequent calls to an instance procedure within a nested 
1253      procedure is not able to access instance variables / procedures. 
1254 
1255      See ambitiousness for a more general workaround. There and elsewhere we
1256      assign *this* (the decision algorithm instance) to *self* and then pass
1257      in self as the first parameter. The receiving instance procedure then
1258      makes sure to call any instance variables and procedures on self rather 
1259      than their implicit and buggy this. The result may be likened to Python's
1260      way of handling instance methods.
1261   */
1262   prob_distance := procedure(p, cp, pp1, pp2) {
1263     /* First pass: (but doesn't handle undefineds and logical equivalents)
1264     return +/ { abs(pp1[v] - pp2[v]) : v in domain(pp1) | v in domain(pp2) };
1265     // Second pass: (but doesn't handle double-counting)
1266     return +/ { abs(pp1[v] - pp2[v2]) : [v, v2] in domain(pp1) >< domain(pp2)
1267                                       | v == v2 || equiv_p_expr(v, v2) };     
1268        Consider three people where the "true" probability is 90%
1269          Person 1: A & B 80%, B & A 70%
1270          Person 2: A & B 80%
1271          Person 3: A & B 80%, B & A 80%
1272        Averaging probabilities in logical equivalents favors person 2 and 3 
1273        over 1 but not 2 over 3.
1274     */
1275     return +/ { abs( pp_equiv(p, cp, pp1)[p_class] -
1276                      pp_equiv(p, cp, pp2)[p_class] )
1277                 : p_class in p_equiv_classes(p, cp)
1278                 | pp_equiv(p, cp, pp1)[p_class] != om &&
1279                   pp_equiv(p, cp, pp2)[p_class] != om };
1280   };
1281 
1282     // Equivalence classes of probability variables based on logical equivalence
1283     p_equiv_classes := procedure(p, cp) {
1284       return { { x : x in range(p) + range(cp)
1285                    | equiv_p_expr(p, cp, p_var, x) }
1286                : p_var in range(p) + range(cp) };
1287     };
1288 
1289     // Converts a possibility probability distribution pp to one that ranges
1290     // over equivalence classes of logically equivalent p vars instead of 
1291     // directly over the p vars.
1292     pp_equiv := procedure(p, cp, pp) {
1293       return { [p_class, avg_p_class(pp, p_class)]
1294                : p_class in p_equiv_classes(p, cp)
1295                | avg_p_class(pp, p_class) != om };
1296     };
1297 
1298       // Return the average of the values pp assigns to the p vars in p_subset
1299       avg_p_class := procedure(pp, p_subset) {
1300         p_vals := [ pp[p_var] : p_var in domain(pp) | p_var in p_subset ];
1301         if (#p_vals == 0) {
1302           return om;
1303         } else {
1304           return avg(p_vals);
1305         }
1306       };
1307 
1308 // Helpers 
1309 // -------
1310 
1311   poss_io_states := cachedProcedure() {
1312     return possible_states(i + o, r);
1313   };
1314 
1315   poss_s_states := cachedProcedure() {
1316     return possible_states(s(), r);
1317   };
1318 
1319   poss_m_states := cachedProcedure() {
1320     return possible_states(m, r);
1321   };
1322 
1323   poss_m2_states := cachedProcedure() {
1324     return possible_states(m2, r);
1325   };
1326 
1327   poss_ms_states := cachedProcedure() {
1328     return possible_states(m + m2, r);
1329   };
1330 
1331   poss_i_states := cachedProcedure() {
1332     return possible_states(i, r);
1333   };
1334 
1335   poss_is_states := cachedProcedure() {
1336     return possible_states(i + s(), r);
1337   };
1338 
1339   poss_o_states := cachedProcedure() {
1340     return possible_states(o, r);
1341   };
1342 
1343   poss_u_states := cachedProcedure() {
1344     return possible_states(range(u), r);
1345   };
1346 
1347 
1348 // Best Compression
1349 // ================
1350 
1351   /* Complexity of Intentional Explanation
1352      -------------------------------------
1353      This can be seen as a more realist and more formal version of Dennett's
1354      Intentional Stance. 
1355        See e.g. https://sites.google.com/site/minddict/intentional-stance-the
1356 
1357      On a first pass, the best intentional explanation can be taken to be the
1358      one which best compresses the brain's behavior.
1359        See better_explanation for further considerations.
1360 
1361      Also see https://en.wikipedia.org/wiki/Minimum_description_length:
1362        "The minimum description length (MDL) principle is a formalization of 
1363         Occam's razor in which the best hypothesis (a model and its parameters) 
1364         for a given set of data is the one that leads to the best compression 
1365         of the data." 
1366 
1367      Since b and d already contain information on their transitions from any 
1368      possible state, we don't need to run this for each time step.
1369   */
1370   complexity := cachedProcedure(b) {
1371     // Complexity of b given d and f + Complexity of d and f
1372     raw_score := k_given(b.t_to_str(), df_to_str(b)) + k(df_to_str(b));
1373     /* Now normalize. Worst case is it doesn't help encode b any better and 
1374        it's just as complex as b. Best case is if it could take 0 bits.  */
1375     worst_score := 2 * k(b.t_to_str());
1376     return 1 - raw_score / worst_score;
1377   };
1378 
1379   to_s := cachedProcedure() {
1380     return str([i, p, cp, u, m, m2, n, o, r, t()]);
1381   };
1382 
1383   f_to_s := procedure(b) {
1384     return str(f(b));
1385   };
1386 
1387   df_to_str := procedure(b) {
1388     return to_s() + "@" + f_to_s(b);
1389     // '@' is just a separator, an otherwise unused character
1390   };
1391 
1392   /* Ambitiousness
1393      -------------
1394      If we only optimized for the best compression of the brain, we face the 
1395      following problem. Suppose the syntax we use is somewhat inefficient. 
1396      Then, the smallest K(b|df) + K(df) might be produced by an empty d with 
1397      all the work being done by a turing machine which generates b using 
1398      essentially the same decision algorithm we would have wanted to use as d 
1399      but with a more efficient syntax. Or it may only put some into the 
1400      official d and sneak in the rest in the turing machine. It seems what
1401      we'd need is something like a not-a-decision-algorithmic-explanation
1402      predicate to apply to the turing machine. Or we want d to be ambitious, 
1403      ideally packing as much of a decision-algorithmic explanation as possible
1404      or at least positively weighing its ambitiousness against potential 
1405      increases in complexity.
1406      
1407      For a given d, we consider the auxiliary hypothesis h, which together w/
1408      d best explains b. We then look for a decomposition of h into [h_i, h_j] 
1409      where [h_i, h_j] is similar to h, and [h_i, h_j] makes d least ambitious, 
1410      meaning that the d* which is a superset of d and maximizes coherence and 
1411      similarity to d + h_i, improves coherence the least with the greatest 
1412      complexity cost. d's ambitiousness in the worst case decomposition is d's 
1413      ambitiousness. In other words, the more there is a d* which is highly 
1414      similar to d + some h_i and very coherent, the less ambitious d is. It 
1415      would leave out too much that could and should be explained using a 
1416      decision algorithmic explanation. If d is very ambitious, there's more 
1417      like just noise left for the auxiliary hypothesis or you'd have to stray 
1418      farther from the auxiliary hypothesis h to subsume it under a decision 
1419      algorithmic explanation.
1420   */
1421   ambitiousness := procedure(b) {
1422     poss_hs := all_strs_lte(#b.t_to_str()) >< all_strs_lte(#b.t_to_str());
1423     self := this; // workaround for setlx bug described in prob_distance
1424     poss_hs := sort_set(poss_hs,
1425       /* Procedure for worse decomposition. Calling it worse rather than
1426          better because we want the [h1, h2] which represents the worst case
1427          scenario as far as ambitiousness goes.  
1428       */
1429       procedure(hs_i, hs_j) {
1430         return hs_score(self, b, hs_i) > hs_score(self, b, hs_j);
1431       });
1432     hs := poss_hs[1];
1433     return ambitiousness_with_hs(b, hs);
1434   };
1435 
1436     hs_score := procedure(self, b, hs_k) {
1437       return self.hs_similarity(b, hs_k) - self.ambitiousness_with_hs(b, hs_k);
1438     };
1439 
1440     hs_similarity := procedure(b, hs_k) { // hs_k == [h1, h2]
1441       return ( k_given(h(b), str(hs_k)) + k_given(str(hs_k), h(b)) ) * -1;
1442       // or ** -1 ?
1443     };
1444 
1445   // Stringified auxiliary hypothesis which together with d explains b.t_to_str
1446   h := procedure(b) {
1447     poss_hs := all_strs_lte(#b.t_to_str());
1448     self := this; // workaround for setlx bug described in prob_distance
1449     poss_hs := sort_set(poss_hs, procedure(h_i, h_j) { // better h explanation
1450                  return h_complexity(self, b, h_i) < h_complexity(self, b, h_j);
1451                });
1452     return poss_hs[1];
1453   };
1454 
1455     h_complexity := procedure(self, b, h_k) {
1456       return k_given(b.t_to_str(), self.df_to_str(b) + "|" + h_k);
1457     };
1458 
1459   ambitiousness_with_hs := procedure(b, hs_k) {
1460     poss_d_stars := { poss_d_star : poss_d_star in supersets_lte(b)
1461                       | poss_d_star.is_valid() &&
1462                         poss_d_star.is_implemented_by(b) };
1463     self := this; // workaround for setlx bug described in prob_distance
1464     poss_d_stars := sort_set(poss_d_stars,
1465       procedure(d_star_i, d_star_j) { // better_d_star
1466         return d_star_sim(self, b, hs_k, d_star_i) - d_star_i.incoherence(b) >
1467                d_star_sim(self, b, hs_k, d_star_j) - d_star_j.incoherence(b);
1468       });
1469     d_star := poss_d_stars[1];
1470 
1471     delta_incoh := ( d_star.incoherence(b) - incoherence(b) ) / incoherence(b);
1472     delta_compl := (k(d_star.df_to_str(b)) - k(df_to_str(b))) / k(df_to_str(b));
1473     /* More incoherent -> more ambitious
1474          because it's like less of real Decision Algorithmic stuff left in h1
1475        Larger complexity -> less ambitious 
1476          because there's more Decision Algorithmic stuff left in hs_1 & not in d
1477          but maybe replace complexity with more like raw size? #states?  */
1478     return delta_incoh - delta_compl;
1479   };
1480 
1481     d_star_sim := procedure(self, b, hs_k, d_star_k) {
1482       raw_sim := (  // roughly, K(d*|d,h1) + K(d,h1|d*_s)
1483         k_given( d_star_k.df_to_str(b), self.df_to_str(b) + "|" + hs_k[1]) +
1484         k_given( self.df_to_str(b) + "|" + hs_k[1],
1485                  d_star_k.simplest_version(b).df_to_str(b) )
1486       );
1487       /* Normalize. Best similarity is if 0 extra bits are required. 
1488          Worst is K(d*) + K(d); the given bits didn't help at all. */
1489       worst_sim := k(d_star_k.df_to_str(b)) + k(self.df_to_str(b));
1490       return 1 - raw_sim / worst_sim;
1491     };
1492 
1493   /* Returns decision algorithms with # of possible states <= #b.poss_states
1494      which are "supersets" of self. It has additional variables and therefore
1495      states, but when considering just the variables shared with self, it is 
1496      isomorphic to self.  */
1497   supersets_lte := procedure(b) {
1498     return { d : d in decision_algorithm.all_lte(#b.poss_states())
1499                | is_subset_of(d) };
1500   };
1501 
1502   // Returns whether this decision algorithm is a "subset" of the given d
1503   is_subset_of := procedure(d) {
1504     // i, o, p, cp, u, m, m2, n, r, t_s, t_o
1505     var_pairs := { [i, d.i], [o, d.o], [p, d.p], [cp, d.cp],
1506                    [u, d.u], [m, d.m], [m2, d.m2] };
1507     vars_are_subsets := forall([vs, d_vs] in var_pairs | vs < d_vs);
1508     if (!vars_are_subsets || #n > #d.n) { return false; }
1509 
1510     n_are_subsets := forall(j in [1..#n] |
1511                        forall(ltr in {'e','m2','o'} |
1512                          n[j][ltr] < d.n[j][ltr] ));
1513     if (!n_are_subsets) { return false; }
1514 
1515     /* For all transitions from a i+s state to a s+o state in self, it is the
1516        case that any transition in d starting from a superset of i+s ends in
1517        a d state that is a superset of the s+o state.  */
1518     return forall([is, so] in t() |
1519              forall(ccs in compatible_complete_states(is, d.poss_is_states()) |
1520                compat(so, d.t()[ccs])
1521              )
1522            );
1523   };
1524 
1525   is_simplest_version := procedure(b) {
1526     return is_equal_to(simplest_version(b));
1527   };
1528 
1529   /* Using simplest_version handles worries with using m as arbitrary data
1530      to smuggle in unlabelled decision algorithmic information. Only the 
1531      decision algorithmic aspects of a decision algorithm should do the 
1532      explaining. This should also handle worries of higher-order utility
1533      information being smuggled in within the memory variables.  */
1534   simplest_version := procedure(b) {
1535     alt_ds := [ d : d in decision_algorithm.all(b)
1536                   | is_isomorphic_to(d) && d.is_implemented_by(b) ];
1537     alt_ds := sort_list(alt_ds, decision_algorithm.simpler);
1538     return alt_ds[1];
1539   };
1540 
1541   is_equal_to := procedure(d) {
1542     l1 := [  i,   o,   p,   cp,   u,   m,   m2,   n,   r,   t_s,   t_o];
1543     l2 := [d.i, d.o, d.p, d.cp, d.u, d.m, d.m2, d.n, d.r, d.t_s, d.t_o];
1544     return l1 == l2;
1545   };
1546 
1547   is_isomorphic_to := procedure(d) {
1548     /* Everything but the memory variables is the same 
1549        And there is an isomorphism between m+m2 states and d.m+d.m2 states
1550        Todo: Also add n minus m2s in n to simple equality comparison.
1551     */
1552     l1 := [  i,   o,   x,   p,   cp,   u,   r,   t_s,   t_o];
1553     l2 := [d.i, d.o, d.x, d.p, d.cp, d.u, d.r, d.t_s, d.t_o];
1554     return l1 == l2 && exists(fm in possible_fms(d) | fm_commutes(fm, d));
1555   };
1556 
1557     // Todo: Also do higher-order m2s
1558     possible_fms := procedure(d) {
1559       /* Construct the space of functions from states of m+m2 to d.m+d.m2,
1560          respecting the m / m2 boundary */
1561       fm_ms  := {
1562         { [ m_st + m2_st,  fm_m[m_st] + fm_m2[m2_st] ]
1563           : [m_st, m2_st] in poss_m_states() >< poss_m2_states() }
1564         : [fm_m, fm_m2] in
1565           function_space(possible_states(m,  r), possible_states(d.m,  d.r)) ><
1566           function_space(possible_states(m2, r), possible_states(d.m2, d.r))
1567       };
1568 
1569       // For convenience map everything else onto its counterpart in d
1570       fms := {
1571         { [iso, apply_fm_m(fm_m, iso) ] : iso in possible_states(i+s()+o, r) }
1572         : fm_m in fm_ms
1573       };
1574       return fms;
1575     };
1576 
1577       // Todo: Also do higher-order m2s.
1578       apply_fm_m := procedure(fm_m, iso) {
1579         iso_m := { [var, val] : [var, val] in iso | var in m || var in m2 };
1580         return { [v, iso[v]] : v in i + s() + o - m - m2 } + fm_m[iso_m];
1581       };
1582 
1583     fm_commutes := procedure(fm, d) {
1584       // Restrict fm to just i + s() variables
1585       fm_is := { [{ [var, val] : [var, val] in iso1 | var in i + s() },
1586                   { [var, val] : [var, val] in iso2 | var in d.i + d.s() }]
1587                     : [iso1, iso2] in fm };
1588       // Restrict fm to just s() + o variables
1589       fm_so := { [{ [var, val] : [var, val] in iso1 | var in s() + o },
1590                   { [var, val] : [var, val] in iso2 | var in d.s() + d.o }]
1591                     : [iso1, iso2] in fm };
1592       return { fm_so[ t()[g_is] ] == t()[ fm_is[g_is] ] : [g_is, h_so] in t()
1593              } == { true };
1594     };
1595 
1596 // Rationality
1597 // ===========
1598 
1599   /* Rational Utility Function
1600      -------------------------
1601      Takes the world w and brain b and collects all possible continuations of
1602      brain/decision states using bs_msr_paths. Then it weights them by 
1603      agential_identity, optimal_rxs, and subjective likelihood and finds the 
1604      center-of-mass / social-choice of the resulting first order utility 
1605      functions. 
1606 
1607      This comes out of the metaethical view developed by Howard Nye and June
1608      Ku and written up at http://www.metaethical.ai/norm_descriptivism.pdf. 
1609      There we focused on the simple case in which the norms we accept or our
1610      higher-order preferences formed a hierarchy with a single fundamental norm 
1611      or highest-order utility function governing our lower-order norms /
1612      functions. (Also see the comments above n := n towards the top.) In that 
1613      case, what we have most reason to do would be what that highest order 
1614      utility function would prescribe, perhaps through the prescription of 
1615      certain lower-level utilty functions.
1616 
1617      Elsewhere, we have acknowledged that things might be more complicated. We 
1618      might have a network of accepted norms / higher-order preferences in which 
1619      they all influence one another rather than a strict top-down hierarchy.
1620      We have suggested that in that case, the weights of each accepted norm and
1621      the parameters governing the network behavior could be treated as the 
1622      fundamental norm. Developing this further, we might model it as a markov 
1623      model with states of norm acceptances influencing other acceptances and 
1624      look for a stationary distribution, i.e. an equilibrium state in which the 
1625      state it transitions to is the same state it started with. 
1626 
1627      However, even that makes a simplifying assumption that the network 
1628      topology of which norms influence which others remains constant. But if 
1629      we imagine a brain's behavior in response to varying inputs, it seems
1630      likely that sometimes, which higher-order preferences influence which 
1631      others would also vary. To handle this, we implement a parliamentary model
1632      to aggregate the various influences across different possible 
1633      continuations of the agent. 
1634        This is inspired by Nick Bostrom and Toby Ord on Parliamentary Models 
1635      for Moral Uncertainty. Although it's applied in a different context, it's
1636      similar in that both are aggregating many utility functions into one.
1637      http://www.overcomingbias.com/2009/01/moral-uncertainty-towards-a-solution.html
1638 
1639      In the end, this does bring us closer to the ideal response / advisor
1640      family of metaethical theories but I think our focus on the influence of
1641      optimal normative judgments remains distinctive and arguably avoids some
1642      worries. We are reading off the first-order utility function from ideal
1643      selves rather than the actions so we have less distortion from the 
1644      idealization process to worry about since there's less to change. We
1645      address remaining worries about distortion with agential identity to 
1646      measure how much one remains the same agent. 
1647 
1648      The idealization is mostly done by weighting more heavily selves who are
1649      the result of normative judgments that best satisfy the decision criteria
1650      governing them. We also incorporate subjective likelihood of the inputs
1651      determining possible continuations to factor in normal operating 
1652      conditions without letting in external factors and violating supervenience
1653      on just the brain.
1654   */
1655   ruf := procedure(w, b) {
1656     t_z := 10 ** 9; // Number of time steps of continuations to consider
1657     /* Todo: Find a more principled way
1658          * "Ask" each branch how much longer to continue?
1659          * Compute infinitely but in a way that asymptotes? 
1660          * Keep computing until it converges, with minimum time? 
1661              If it never converges? */
1662     voter_wts := {};
1663     t_i := 1;
1664     bs1 := drop_time(b.a[b.n]);
1665     while (t_i <= t_z) {
1666       bmps := bs_msr_paths(w, b, bs1, t_i);
1667       voter_wts := voter_wts + {
1668         [ additive_to_ordinal(f(b)[first(last(bmp))], 0), // voter util func
1669           time_wt(t_i) * voter_weight(w, b, bmp) ]        // voter weight
1670         : bmp in last(bmps)
1671       };
1672       t_i := t_i + 1;
1673     }
1674     // Poss social choice utility functions
1675     psc_u_dom := domain(p) + domain(u);
1676     psc_u_rng := r[first(range(u))] + { om };
1677     psc_us := choices({ { [expr, r_i] : r_i in psc_u_rng }
1678                         : expr in psc_u_dom });
1679     psc_us := { { [expr, r_i] : [expr, r_i] in psc_u | r_i != om }
1680                 : psc_u in psc_us }; // filter out om
1681     return sort_set(psc_us,
1682                     procedure(psc_u1, psc_u2) {
1683                       return rat_psc_dist(voter_wts, psc_u1) <
1684                              rat_psc_dist(voter_wts, psc_u2);
1685                     })[1];
1686   };
1687 
1688     time_wt := procedure(t_i) { return 1; }; // Placeholder constant function
1689 
1690     voter_weight := procedure(w, b, bmp) {
1691       msr := last(last(bmp));
1692       bss := [ bs : [bs, msr] in bmp ]; // Brain StateS
1693 
1694       cached_ref := om;
1695       wb := new_wb(w, b, bss);
1696       w2 := first(wb);
1697       b2 := last(wb);
1698 
1699       dss := [ f(b2)[bs] : bs in bss ];
1700       agential_id := agential_identity(dss, w2, b2);
1701       return msr * agential_id * optimal_rxs(w2, b2, dss);
1702     };
1703 
1704     // Possible Social Choice for Rationality (as opposed to in optimal_rxs)
1705     rat_psc_dist := procedure(voter_wts, psc_u) {
1706       ord_u := card_u_to_ord_u( add_u_to_card_u(psc_u) );
1707       return +/ { voter_wts[voter] * ord_u_dist(voter, ord_u) ** 2
1708                   : voter in domain(voter_wts) };
1709     };
1710 
1711   new_wb := procedure(w, b, bss) {
1712     b2   := b;
1713     b2.a := b2.a + [assign_time(bss[t_j], b.n + t_j) : t_j in [1..#bss]];
1714     b2.n := #(b2.a);
1715 
1716     w2   := w;
1717     w2.a := w2.a + [assign_time(bss[t_j], b.n + t_j) : t_j in [1..#bss]];
1718     w2.n := #(w2.a);
1719 
1720     return [w, b];
1721   };
1722 
1723   /* Optimal Normative Judgments / Prescriptions
1724      -------------------------------------------
1725      Returns a score for how much the list of decision states has had outputs
1726      of higher-order utility functions (i.e. normative judgments / 
1727      prescriptions) that maximized higher-order utility. If the state space
1728      of prescriptions (n[j]['o']) overlap, we favor those states voted for by 
1729      more accepted norms through a synchronic social choice function.
1730   */
1731   optimal_rxs := procedure(w, b, dss) {
1732     opts := [ optimality(w, b, dss[j], dss[j+1], b.n - (#dss - j) )
1733               : j in [1..(#dss - 1)] ];
1734     return 0.5 * min(opts) + 0.5 * avg(opts);
1735   };
1736 
1737     optimality := procedure(w, b, ds, ds2, t_i) {
1738       ord_us := { additive_to_ordinal(ds, j) : j in [1..#n] };
1739 
1740       /* Social choice to balance multiple synchronic utility functions if
1741          they conflict.  */
1742       psc_base_exprs := +/ { ue_base_exprs(j) : j in [1..#n] };
1743       psc_state_space := choices({ {[subform, true], [subform, false]}
1744                                    : subform in psc_base_exprs });
1745       sync_psc_us := choices({ { [state, r_i] : r_i in r[first(range(u))] }
1746                                : state in psc_state_space });
1747       sc_u := sort_set(psc_us,
1748                        procedure(psc_u1, psc_u2) {
1749                          return sync_psc_dist(ord_us, psc_u1) <
1750                                 sync_psc_dist(ord_us, psc_u2);
1751                        })[1];
1752 
1753       n_o_vars := +/ { n_i['o'] : n_i in n }; // Prescription vars
1754       // Observed state of those vars in ds2
1755       obs_rxs := { [n_o_var, ds2[n_o_var]] : n_o_var in n_o_vars };
1756 
1757       // Collect all possible prescriptions along with their (higher) utilities.
1758       util_rxs := {
1759         [rxs, util_of_rxs(w, b, t_i, ds, ds2, psc_base_exprs, sc_u, rxs)]
1760         : rxs in poss_rx_states()
1761       };
1762       // Sort them by optimality.
1763       util_rxs := sort_set(util_rxs, procedure(rxs_u1, rxs_u2) {
1764                                        return last(rxs_u1) > last(rxs_u2);
1765                                      });
1766       // Return a normalized score for the obs_rxs by their percentile rank.
1767       return 1 - index_of(
1768         [ obs_rxs,
1769           util_of_rxs(w, b, t_i, ds, ds2, psc_base_exprs, sc_u, obs_rxs) ],
1770         util_rxs
1771       ) / #util_rxs;
1772     }; // end optimality
1773 
1774       sync_psc_dist := procedure(ord_us, sync_psc_u) {
1775         ord_psc_u := card_u_to_ord_u(sync_psc_u);
1776         return +/ { ord_u_dist(ord_u, ord_psc_u) ** 2 : ord_u in ord_us };
1777       };
1778 
1779       util_of_rxs := procedure(w, b, t_i, ds, ds2, psc_base_exprs, sc_u, rxs) {
1780         /* For the sake of simplicity, we assume that higher-order utilities 
1781            depend on only intrinsic properties of brains / decision algorithms.
1782            This means that being given the future brain states bss is 
1783            sufficient to tell us the utility of various prescriptions under 
1784            this assumed future. This assumption can be quite plausible for 
1785            higher-order utilities that enforce consistency, coherence or some
1786            form of reflective equilibrium. If we abandon this assumption, we
1787            may need to model a set of possible worlds and a probability 
1788            distribution or measure over them.
1789 
1790            We might typically expect higher-order utilities to be placed in 
1791            some change in lower-order utility and then a recognition that some
1792            prescription would more-or-less immediately cause that change would 
1793            cause the prescription and in turn, the intended change. But since
1794            we can neither assume that the higher-order utility is placed in 
1795            a specific, concrete lower-order utility change in the near future as
1796            opposed to more open-ended properties nor assume the immediacy of the
1797            prescription's causation of a lower-order utility change, we just 
1798            look very far into the future and sum up the utility of any events 
1799            in that 4-dimensional space-time block that would be included in the
1800            reference of higher-order utility expressions. Since the network
1801            topology of the higher-order utility functions may be path dependent,
1802            we also use the bs_msr_paths to aggregate over the different futures
1803            made possible by various input states weighted by their subjective
1804            likelihood.
1805         */
1806         new_ds2 := rxs + { [var, val] : [var, val] in ds2
1807                                       | !(var in domain(rxs)) };
1808         bs := inv(f(b))[new_ds2];
1809         bmps := bs_msr_paths(w, b, bs, 10 ** 9);
1810         future_wts := {                      // collect a ... 
1811           [ [first(bs_msr) : bs_msr in bmp], // list of brain states, bss,
1812             last(last(bmp)) ]                // & measure from the final state
1813           : bmp in last(bmps)                // for each bs_msr_path
1814         }; // to get { [bss, msr], ... }
1815 
1816         return +/ {
1817           msr * util_of_future(w, b, t_i, ds, psc_base_exprs, sc_u, bss)
1818           : [bss, msr] in future_wts
1819         };
1820       }; // end util_of_rxs
1821 
1822         util_of_future := procedure(w, b, t_i, ds, psc_base_exprs, sc_u, bss) {
1823           wb := new_wb(w, b, bss);
1824           w2 := first(wb);
1825           b2 := last(wb);
1826 
1827           future_state := { [expr, eval_expr(w2, b2, t_i, ds, bss, expr)]
1828                             : expr in psc_base_exprs };
1829           return sc_u[future_state];
1830         };
1831 
1832         // Returns a boolean for whether the expr is true in that future
1833         eval_expr := procedure(w, b, t_i, ds, bss, expr) {
1834           p_var := p[expr];
1835           p_event := { [p_var, ds[p_var]] };
1836           b_event := first( lf_inv_at(b,t_i)[p_event] );
1837           ref_b_event := ref(w,b)[b_event];
1838           disjuncts := { restricted_eq(w.actual_state(), disjunct)
1839                          : disjunct in ref_b_event };
1840           return #{ disj : disj in disjuncts | disj == true } > 0;
1841         };
1842 
1843   /* Extensional rational utility */
1844   ext_ruf := procedure(w, b) {
1845     uf := ruf(w,b);
1846     // Map expressions to their referents.
1847     expr_ref := { [expr, ref_expr(w, b, b.n, expr)] : [expr, r_i] in uf };
1848     // Map every possible world state to its utility.
1849     return { [w_state, w_util(uf, expr_ref, w_state)]
1850              : w_state in w.cm().poss_states() };
1851   };
1852 
1853     w_util := procedure(uf, expr_ref, w_state) {
1854       /* Truths, or utility relevant truths, I suppose. Collect the expressions 
1855          in the rational utility function that had this w_state among its set of
1856          truthmakers / referents.  */
1857       truths := { expr : expr in domain(expr_ref) | w_state in expr_ref[expr] };
1858       // Sum the utilities 
1859       return +/ { uf[expr] : expr in truths };
1860     };
1861 
1862   /* Agential Identity
1863      -----------------
1864      Agential Identity is much like Personal Identity a la Derek Parfit's 
1865      Reasons and Persons. But where personal identity roughly concerns when 
1866      someone is the same person for the purposes of care, agential identity
1867      concerns when someone is the same agent for purposes of trusting their
1868      normative judgments. 
1869 
1870      If we are thinking in terms of ideal advisor / observer theories, we might
1871      be considering hypothetical versions of ourselves with more information 
1872      or other qualities. But consider the full knowledge of what heroin feels 
1873      like. If we had that information, we might become addicts and change our 
1874      values. But this probably isn't a change that would make us want to defer
1875      to that self's values. We may protest that something in the process of
1876      adding this knowledge has changed our values, not illuminated them.
1877      
1878      The concept of agential identity is meant to measure how much other 
1879      possible selves truly reflect *our* values and agency. Generally, this 
1880      will be measured either by similarity to our current selves or by
1881      identifying the changes as cognitive changes made through inference or
1882      applying higher-order decision criteria to values. 
1883   */
1884   agential_identity := procedure(dss, w, b) { // dss : list of decision states
1885     return 0.5 * agential_continuity(dss, w, b) +
1886            // kinda like: 0.5 * agential_connection(dss[1], {}, dss[#dss])
1887            0.5 * (0.5 * p_conn(dss[1], last(dss)) +
1888                   0.5 * u_conn(dss[1], {}, last(#dss)));
1889   };
1890   agential_continuity := procedure(dss, w, b) {
1891     conns := [agential_connection(dss[j], dss[j+1], dss[j+2], w, b, j)
1892               : j in [1..(#dss - 2)] ];
1893     return 0.5 * min(conns) + 0.5 * avg(conns);
1894   };
1895   agential_connection := procedure(ds1, ds2, ds3, w, b, t_i) {
1896     return 0.5 * avg([p_conn(ds1, ds2), p_conn(ds2, ds3)]) +
1897            0.5 * u_conn(ds1, ds2, ds3, w, b, t_i);
1898   };
1899   /* Agential continuity restricted to p_conn.
1900      Used for measuring diachronic rationality in incoherence. */
1901   p_conn_continuity := procedure(dss) {
1902     return avg({ p_conn(dss[j], dss[j+1]) : j in [1..(#dss - 1)] });
1903   };
1904   /* Measures connectedness of credences between two decision states. They
1905      are connected by either having similar probability values or by reflecting
1906      cognitive changes, i.e. updating using conditional probabilities and
1907      resulting in probability values that are close to (locally) ideally 
1908      rational inferences. */
1909   p_conn := procedure(ds1, ds2) {
1910     return avg({
1911       // Gain full credit of the higher of similarity or cognitive score and
1912       max( sim_cog_scores(ds1, ds2, p_i) ) +
1913         // Close the remaining distance proportional to
1914         (1 - max( sim_cog_scores(ds1, ds2, p_i) )) *
1915         // the other of the similarity or cognitive score
1916         min( sim_cog_scores(ds1, ds2, p_i) )
1917       // For every p var
1918       : p_i in p_vars()
1919     });
1920   };
1921 
1922     sim_cog_scores := procedure(ds1, ds2, p_i) {
1923       return [p_i_sim(ds1, ds2, p_i), cog_change(ds1, ds2, p_i)];
1924     };
1925 
1926     // Similarity / closeness. Returns 0-1.
1927     p_i_sim := procedure(ds1, ds2, p_i) {
1928       max_dist := max({ abs(ds1[p_i] - r_i) : r_i in r[p_i] });
1929       return 1 - abs(ds1[p_i] - ds2[id_var(p_i)]) / max_dist;
1930     };
1931 
1932     // Cognitive change. Returns 0-1.
1933     cog_change := procedure(ds1, ds2, p_i) {
1934       if (p_i in range(p)) { // rather than range(cp)
1935         return 0;
1936       } else {
1937         arr := p_inv(p_i);
1938         hyp := arr[1];
1939         evd := arr[2];
1940         ancs := non_shared_ancestors(p_i, m);
1941         if (p[hyp] in ancs && p[evd] in ancs && cp[[evd, hyp]] in ancs) {
1942           /* For simplicity we just use the p vals in the immediate parent 
1943              (ds1). Really we should construct and use 
1944              non_shared_ancestors_with_depth and pull the p values from the 
1945              correct time slice. Similarly for id_expr. */
1946 
1947           // Rational P(H|E) = P(E|H) * P(H) / P(E)
1948           rat_h_e := ds1[cp[[evd, hyp]]] * ds1[p[hyp]] / ds1[p[evd]];
1949           // Distance from rational
1950           dist_rat := abs(rat_h_e - ds2[ cp[ [id_expr(hyp), id_expr(evd)] ] ]);
1951           // Maximum distance possible
1952           max_dist := max({ abs(rat_h_e - r_i) : r_i in r[p_i] });
1953           return 1 - dist_rat / max_dist;
1954         } else {
1955           return 0;
1956         }
1957       } // end else 
1958     }; // end cog_change
1959 
1960   /* Takes an expression e and returns the expression at the next time step
1961      that would be identical. In many cases, this is just the same expression e,
1962      but we define input variables as indexicals so the same expression next
1963      time step does not refer to the expression in this time step. So, we
1964      introduce a new term At(expr, -1) which at that time step, would refer to 
1965      the same as this time's expr.  
1966   */
1967   id_expr := procedure(expr) {
1968     if (isList(expr)) { // in case it's called on p_inv of cp var
1969       return [ id_expr(expr_i) : expr_i in expr ];
1970     } else {
1971       if (p[expr] in i) { // current input
1972         return At_t(expr, -1);
1973       } else if (isTerm(expr) && fct(expr) == "At_t") { // past, add 1 more step
1974         return At_t(args(expr)[1], args(expr)[2] - 1);
1975       } else {
1976         return expr;
1977       }
1978     }
1979   };
1980 
1981   /* Takes a prob variable and finds its expr preimage, then finds the expr that
1982      in the next time would be identical to this expr. Returns that expr var. */
1983   id_var := procedure(p_var) {
1984     return p[id_expr(p_inv(p_var))];
1985   };
1986 
1987   poss_i_events := procedure() {
1988     return +/ { possible_states(isub, r) : isub in ne_pow(i) };
1989   };
1990 
1991   e_vars := procedure() {
1992     return +/ { range(n_i['e']) : n_i in n };
1993   };
1994 
1995   poss_ue_events := procedure() {
1996     return +/ { possible_states(ue_sub, r)
1997                 : ue_sub in ne_pow(range(u) + e_vars()) };
1998   };
1999 
2000   poss_p_i_events := procedure() {
2001     return +/ { { [p_i, r_i] : r_i in r[p_i] } : p_i in range(p) };
2002   };
2003 
2004   poss_p_var_events := procedure() {
2005     return +/ { { [p_i, r_i] : r_i in r[p_i] } : p_i in p_vars() };
2006   };
2007 
2008   // Possible states of all prescription vars
2009   poss_rx_states := procedure() {
2010     n_o_vars := +/ { n_i['o'] : n_i in n };
2011     return possible_states(n_o_vars, r);
2012   };
2013 
2014   // Possible single prescriptions for some higher-order utility
2015   poss_rxs := procedure() {
2016     n_os := { n_i['o'] : n_i in n };
2017     return +/ { possible_states(n_o, r) : n_o in n_os };
2018   };
2019 
2020   /* Utility Function Connectedness
2021      ------------------------------
2022      Measures how much the utility function has either remained similar or 
2023      changed as the result of a cognitive process involving i) recognition that 
2024      some prescription (rx) would cause a utility change, ii) that recognition
2025      in fact causes the prescription and iii) that prescription causes the 
2026      utility change. These are to be contrasted with changes to utilities from
2027      mere noise, associations or biases.  */
2028   u_conn := procedure(w, b, t_i, ds1, ds2, ds3) {
2029     actual_rxs := { rx : rx in poss_rxs()
2030                        | rx_to_tuple(w, b, t_i, ds1, ds2, ds3, rx) != {} };
2031     rx_subsets := pow(actual_rxs);
2032 
2033     dists := { dist_to(w, b, t_i, ds1, ds2, ds3, rx_subset)
2034                : rx_subset in rx_subsets };
2035     return 1 - min(dists);
2036   }; // end u_conn
2037 
2038     /* Takes a prescription and returns the tuples if any which realize i, ii
2039        and iii above.  */
2040     rx_to_tuple := procedure(w, b, t_i, ds1, ds2, ds3, rx) {
2041       poss_tuples := poss_p_var_events() >< poss_i_events()  ><
2042                                p_range() >< poss_ue_events();
2043       return { tuple : tuple in poss_tuples
2044                      | observed(w, b, t_i, rx, ds1, ds2, ds3, tuple) };
2045     };
2046 
2047       observed := procedure(w, b, t_i, rx, ds1, ds2, ds3, tuple) {
2048         f_p_i  := tuple[1]; // Assignment of value to p_i, rx []-> ue
2049         g_i_rx := tuple[2]; /* Assignment of value to a prescription which is
2050                                introspectively accessible as i vars */
2051         p_ue_val := tuple[3]; // Credence placed in ue (the consequent of []->)
2052         h_ue   := tuple[4]; // ue event
2053         return f_p_i != {} && g_i_rx != {} && h_ue != {} &&
2054           // g_i_rx and rx are implemented by the same brain event and 
2055           // the decision events occur at the right times
2056           lf_inv(g_i_rx) == lf_inv(rx) && restricted_eq(ds2, g_i_rx) &&
2057           restricted_eq(ds2, rx) && restricted_eq(ds1, f_p_i) &&
2058           restricted_eq(ds3, h_ue) &&
2059           // first(domain(f_p_i)) is p_i and p_i is _ []-> _
2060           fct(p_var_inv(first(domain(f_p_i)))) == "BoxArrow" &&
2061           /* The ref of the brain event of ds1 implementing the placement of 
2062              prob in the p_i expr == lf_inv h_ue in two time steps
2063              i.e. the consequent of []-> refers to the intended ue change */
2064           /* Todo: Allow for changes that had i, ii, and iii but only partially 
2065              updated the utility / eval functions.  */
2066           ref(w,b)[lf_inv_at(b, t_i)[
2067             {[ p[args(p_var_inv(first(domain(f_p_i))))[2]], p_ue_val ]}
2068           ]] == lf_inv_at(b, t_i + 2)[h_ue] &&
2069           // p_i influences each var in rx
2070           first(domain(f_p_i)) in { parents(v) : v in domain(rx) } &&
2071           // each var in rx influences each var in the utility change 
2072           forall(v in g_i_rx | v in { parents(ue) : ue in domain(h_ue) });
2073       };
2074 
2075     /* Given a subset of observed prescriptions, we compare the resulting 
2076        utility function to what was prescribed within that subset. The overall 
2077        u_conn score will take the minimum distance seen this way.
2078          Todo: Allow for changes that had i, ii, and iii but only partially 
2079        updated the utility / eval functions.  */
2080     dist_to := procedure(w, b, t_i, ds1, ds2, ds3, rx_subset) {
2081       u1 := { [u_i, ds1[u_i]] : u_i in u };
2082       e1 := { u1 } + { { [e_i, ds1[e_i]] : e_i in n_i['e'] } : n_i in n };
2083 
2084       rxd := { [e_i, rx_e_i_val(w, b, t_i, ds1, ds2, ds3, rx_subset, e_i)]
2085                : e_i in domain(e1) };
2086       rxd_u := { [u_i, rxd[u_i]] : u_i in u };
2087       rxd_n := [ { [e_i, rxd[e_i]] : e_i in n_i['e'] } : n_i in n ];
2088       rxd_e := [ rxd_u ] + rxd_n;
2089       // The prescribed util / eval change
2090       ord_rxd_e := [ additive_to_ordinal(rxd_e[j], j) : j in [1..#rxd_e] ];
2091       // The actual change in ds3
2092       ord_e := [ additive_to_ordinal(e_func(ds3, j), j) : j in [1..#rxd_e] ];
2093 
2094       e_dists := [ ord_u_dist(ord_e[j], ord_rxd_e[j]) : j in [1..#ord_e] ];
2095 
2096       // Todo: Add weights? 
2097       return avg(e_dists);
2098     }; // end dist_to
2099 
2100       /* Given a subset of observed prescriptions and an e_i var, return the
2101          value prescribed for the e_i var.  */
2102       rx_e_i_val := procedure(w, b, t_i, ds1, ds2, ds3, rx_subset, e_i) {
2103         /* The prescribed values (rxds) are the values of the e_i var in the 
2104            ue event from rx_to_tuple (for any rxs that prescribe for e_i). */
2105         rxds := { rx_to_tuple(w, b, t_i, ds1, ds2, ds3, rx_i)[4][e_i]
2106           : rx_i in rx_subset  // h_ue[e_i]
2107           | e_i in domain(rx_to_tuple(w, b, t_i, ds1, ds2, ds3, rx_i)[4]) };
2108             // e_i in domain(h_ue)
2109         if (#rxds > 0) {
2110           // If #rxds > 1, they should all be the same anyway bc it was observed
2111           return first(rxds);
2112           /* Todo: In later versions, we probably want to allow for non-exact 
2113              changes and then choose the one closest to the observed.  */
2114         } else {
2115           return e1[e_i]; // No change
2116         }
2117       };
2118 
2119   // Utility / Evaluation Function implemented in a given decision state ds
2120   e_func := procedure(ds, j) {
2121     if (j == 0) {
2122       return { [u_i, ds[u_i]] :  u_i in range(u) };
2123     } else {
2124       return { [e_i, ds[e_i]] : e_i in range(n[j]['e']) };
2125     }
2126   };
2127 
2128   ue := procedure(j) {
2129     if (j == 0) {
2130       return u;
2131     } else {
2132       return n[j]['e'];
2133     }
2134   };
2135 
2136   ue_base_exprs := procedure(j) {
2137     return +/ { subformulas(x_i) : x_i in domain(ue(j)) };
2138   };
2139 
2140   /* Takes a decision state, i.e. a function from (among other things) utility
2141      variables to numbers (utilities), and an index j and returns a set of 
2142      SuccEq (succeeds or equals) setlx functors. If j is 0, then u is converted
2143      to ordinal, otherwise, n[j]['e'] is. */
2144   additive_to_ordinal := procedure(ds, j) {
2145     return add_u_to_card_u(e_func(ds,j));
2146   };
2147 
2148   /* Brain-State Measure Paths 
2149      -------------------------
2150      Returns a list where for each future time, there is a set of 
2151      bs_msr_paths, each of which in turn is a list of bs_msrs for each time 
2152      from the first to the given time. A bs_msr is a pair [bs, msr]. 
2153     
2154      One can think in terms of a tree. At the root is the given starting brain 
2155      state. Each branch is a continuation with different inputs and a measure 
2156      for that branch based on subjective likelihood. A path is a full branch 
2157      line from root to leaf collecting pairs of brain states and measures along
2158      the way. 
2159 
2160      e.g. [bs_msr_paths_t1, ..., bs_msr_paths_tz]
2161        bs_msr_paths_tj := { bs_msrs_tj_1, bs_msr_tj_2, ... }
2162          bs_msrs_tj_k := [ bs_msr_t1, ..., bs_msr_tj ]
2163            bs_msr_tj := [bs_tj, msr]
2164 
2165                      [ bs1, 1 ]
2166                     /   |      \
2167                is1 /    |is2    \  is3   (input states leading to continuations)
2168                   /     |        \  
2169                  /      |         \
2170      [bs_t2a, 1/6]  [bs_t2b, 1/2]  [bs_t2c, 1/3]  
2171       /    |    \     /   |   \      /   |   \
2172      ...
2173 
2174      There's some duplication and space inefficiency in having the previous 
2175      times' brain states and measures present both in the earlier elements of
2176      the top level list as well as within the set of paths for each element.
2177      We probably could do without it but it makes it a little more convenient
2178      when calculating the voter weights for each time's voters in computing
2179      the rational utility function. Some aspects of those weights like the 
2180      agential identity depend on the entire path taken and not just the final
2181      brain state.
2182   */
2183   bs_msr_paths := procedure(w, b, bs1, t_j) {
2184     if (t_j == 1) {
2185       return [{[[bs1, 1]]}];
2186     } else {
2187       return bs_msr_paths(w, b, bs1, t_j - 1) +
2188         [+/ {
2189            { bmp + [bs_msr]
2190              : bs_msr in next_bs_msrs(w, b, last(bmp)[1], last(bmp)[2], t_j-1) }
2191            : bmp in last(bs_msr_paths(w, b, bs1, t_j - 1))
2192          }];
2193     }
2194   };
2195 
2196   /* Takes a brain, brain state, measure and prev t_i and returns the possible
2197      continuations in the next time step together with their measure 
2198      roughly proportional to the subjective likelihood.  */
2199   next_bs_msrs := procedure(w, b, bs_i, msr, t_j) {
2200     bs_j := b.response(bs_i);
2201     // Possible brain input states from each di possible decision input states.
2202     pbis := { [di, lf_inv(b)[di]] : di in poss_i_states() };
2203     /* Possibly empty brain input states since we now filter by compatibility
2204        with the brain response bs_j. (Recall that input vars may overlap with
2205        other vars that bs_j may dictate values for.) */
2206     pe_bis := { [di, { pbi : pbi in pbi_s | compat(pbi, bs_j) }]
2207                 : [di, pbi_s] in pbis };
2208     // pbi_s = possible brain input sets
2209     bis := { [di, pbi_s] : [di, pbi_s] in pe_bis | pbi_s != {} };
2210 
2211     /* For each brain input state of each decision input state, map the total
2212        brain state including the response bs_j to its new measure. 
2213 
2214        Ideally, the new measure would be proportional to the brain input 
2215        state's subjective probability. This was under development but doing it 
2216        in a way that is biologically plausible presents a challenge. To avoid 
2217        an implausible number of explicit probabilities, we'd probably have to 
2218        look at the closest coherent extensions of the current probability 
2219        distribution as well as allow for indirect reference, which introduces 
2220        its own challenges in translating probabilities from direct to indirect.
2221        
2222        Although I by no means think these challenges insurmountable, the 
2223        technical complications involved seemed too great for a first version
2224        and not worth delaying the initial release over. As a placeholder here,
2225        we simply assume each possible decision input state is equally probable
2226        and each brain state implementing a given decision input state is
2227        similarly equally probable to other such implementations.
2228 
2229        While incorporating the subjective probability seems preferable, it's
2230        also worth mentioning that it may not actually be necessary. The main
2231        motivation is to make continuations of the agent under very weird 
2232        scenarios matter less. But if the worry there is that those 
2233        circumstances are distorting the agent's values, we already have the 
2234        agential_identity measure to discount them. In other words, agential
2235        identity may screen off a circumstance's weirdness from its relevance.
2236     */
2237     return +/ {
2238       { [ bi + bs_j, msr * (1 / #domain(bis)) * (1 / #bis[di]) ]
2239         : bi in bis[di] }
2240       : di in domain(bis)
2241     };
2242   };
2243 
2244 
2245 // Class Methods
2246 // =============
2247 
2248   static {
2249     /* Function from a causal_model of a brain b and a set of decision 
2250        algorithm candidates to the decision algorithm that is the best 
2251        compression of the brain.  */
2252     implemented_by := cachedProcedure(b, ds) {
2253       if (ds == {}) {
2254         ds := decision_algorithm.all(b);
2255       }
2256       /* Requiring simplest version helps ensure higher-order ambitiousness,
2257          i.e. prevents smuggling in higher-order norms within the n[j]['m2'] 
2258          variables. In future iterations, this might be softened.  */
2259       ds := [d : d in ds | d.is_implemented_by(b) && d.is_simplest_version(b)];
2260       ds := sort_list(ds, decision_algorithm.better_explanation(b));
2261       return ds[1];
2262     };
2263 
2264     /* Given a causal model of a brain, return the set of possible decision 
2265        algorithms which could be implemented by it based simply on number of 
2266        states. A decision algorithm can't be more fine-grained than its 
2267        substrate.  */
2268     all := cachedProcedure(b) {
2269       return { d : d in decision_algorithm.all_lte(#b.poss_states()) };
2270     };
2271 
2272     // All with #poss_states <= z
2273     // Todo: Add possible At_t exprs and predicates and quantifiers
2274     all_lte := cachedProcedure(z) {
2275       poss_i := da().ltr_up_to('I', z/2);
2276       var_ltrs := { 'A', 'B', 'C' };
2277       poss_base_vars := da().vars_up_to(var_ltrs, n);
2278       poss_exprs := { da().add_up_to(pbe, z/2)
2279                       : pbe in poss_base_vars };
2280       poss_p := da().poss_expr_to_ltr(poss_exprs, 'P');
2281       // e.g. { [[x, y], 'P(x|y)'], ... } 
2282       poss_cp := pow({ [[x, y], 'P(' + str(x) + '|' + str(y) + ')']
2283                        : [x, y] in poss_exprs >< poss_exprs });
2284       poss_o := da().ltr_up_to('O', z/2);
2285       poss_u := da().poss_expr_to_ltr(poss_exprs, 'U');
2286       poss_m := da().ltr_up_to('M', z/2);
2287       poss_m2 := da().ltr_up_to('M2', z/2);
2288       poss_n := +/ { choices({ da().poss_nj(poss_exprs, z, j) : j in [1..x] })
2289                      : x in [1..z] };
2290       return +/ {
2291         +/ { da().poss_ds(d_i, d_o, d_p, d_cp, d_u, d_m, d_m2, d_n, poss_r)
2292              : poss_r in da().poss_rs(z, d_i, d_o, d_p, d_cp,
2293                                          d_u, d_m, d_m2, d_n) }
2294         : [d_i, d_o, d_p, d_cp, d_u, d_m, d_m2, d_n] in
2295           da().cart_prod(poss_i, poss_o, poss_p, poss_cp,
2296                          poss_u, poss_m, poss_m2, poss_n)
2297       };
2298     }; // end all_lte
2299 
2300       ltr_up_to := procedure(ltr, j) {
2301         return { ltr + '_' + y : y in [1..j] };
2302       };
2303 
2304       vars_up_to := procedure(ltrs, j) {
2305         return +/ { da().ltr_up_to(ltr, j) : ltr in ltrs };
2306       };
2307 
2308       add_one := procedure(exprs) {
2309         unary := { 'Not' };
2310         binary := { 'And', 'Or', 'Implies', 'BoxArrow' };
2311         add_un := { makeTerm(op, [expr]) : [op, expr] in unary >< exprs };
2312         add_bi := { makeTerm(op, [expr1, expr2])
2313                     : [op, [expr1, expr2]] in binary >< (exprs >< exprs) };
2314         return exprs + add_un + add_bi;
2315       };
2316 
2317       add_up_to := procedure(base_vars, j) {
2318         if (j == 0) {
2319           return base_vars;
2320         } else {
2321           return da().add_one( da().add_up_to(base_vars, j - 1) );
2322         }
2323       };
2324 
2325       expr_to_ltr_var := procedure(expr, ltr) {
2326         return ltr + '(' + replace(str(expr), '"', '') + ')';
2327       };
2328 
2329       poss_expr_to_ltr := procedure(poss_exprs, ltr) {
2330         return { { [expr, da().expr_to_ltr_var(expr, ltr)] : expr in subset }
2331                  : subset in pow(poss_exprs) };
2332       };
2333 
2334       poss_nj := procedure(poss_exprs, z, j) {
2335         poss_e_nj := da().poss_expr_to_ltr(poss_exprs, 'E_n' + j);
2336         poss_m2_nj := da().ltr_up_to('M2_n' + j, z/2);
2337         poss_o_nj := da().ltr_up_to('O_n' + j, z/2);
2338         return { { ['e', e_n], ['m', m_n], ['o', o_n] }
2339           : [e_n, m_n, o_n] in cart_prod([poss_e_nj, poss_m2_nj, poss_o_nj]) };
2340       };
2341 
2342       poss_rs := procedure(z, d_i, d_o, d_p, d_cp, d_u, d_m, d_m2, d_n) {
2343         poss_p_vals := { { j / x : j in [0..x] } : x in [1..z] };
2344         poss_ue_vals := ne_pow({ j : j in [1..2**z] });
2345         poss_z_vals := { { j : j in [1..x] } : x in [1..z] };
2346         poss_o_vals := poss_z_vals;
2347 
2348         poss_r_ps := {
2349           { [i_var,  p_vals] : i_var  in d_i         } +
2350           { [p_var,  p_vals] : p_var  in range(d_p)  } +
2351           { [cp_var, p_vals] : cp_var in range(d_cp) }
2352           : p_vals in poss_p_vals
2353         };
2354 
2355         return +/ choices(
2356           { poss_r_ps, da().poss_r_ue(poss_ue_vals, d_u) } +
2357           { da().poss_r_ue(poss_ue_vals, n_i['e']) : n_i in d_n } +
2358           { da().poss_r_z(poss_z_vals, d_z) : d_z in { d_o, d_m, d_m2 } }
2359         );
2360       }; // end poss_rs
2361 
2362         poss_r_ue := procedure(poss_ue_vals, ue) {
2363           return { { [ue_var, ue_vals] : ue_var in range(ue) }
2364                    : ue_vals in poss_ue_vals };
2365         };
2366 
2367         poss_r_z := procedure(poss_z_vals, d_z) {
2368           return choices({ {z_var} >< poss_z_vals : z_var in d_z });
2369         };
2370 
2371       poss_ds := procedure(z, d_i, d_o, d_p, d_cp, d_u, d_m, d_m2, d_n, d_r) {
2372         poss_is_state := possible_states(d_i + d_p + d_cp +
2373                                          d_u + d_m + d_m2 + d_n, d_r);
2374         poss_s_state  := possible_states(d_p + d_cp + d_u +
2375                                          d_m + d_m2 + d_n, d_r);
2376         poss_t_s := function_space(poss_is_state, poss_s_state);
2377 
2378         poss_p_u_m2 := possible_states(d_p + d_u + d_m2, d_r);
2379         poss_o_state := possible_states(d_o, d_r);
2380         poss_t_o := function_space(poss_p_u_m2, poss_o_state);
2381 
2382         das := {
2383           decision_algorithm(d_i, d_o, d_p, d_cp, d_u, d_m, d_m2, d_n,
2384                              d_r, d_t_s, d_t_o, om, om, om)
2385           : [d_t_s, d_t_o] in poss_t_s >< poss_t_o
2386         };
2387         return { da : da in das | da.is_valid() && #da.poss_states() <= z};
2388       };
2389 
2390     better_explanation := cachedProcedure(b) {
2391       return procedure(d1, d2) {
2392                return da().expl_score(b, d1) > da().expl_score(b, d2);
2393              };
2394     };
2395 
2396       expl_score := procedure(b, d_i) {
2397         return d_i.ambitiousness(b) - d_i.complexity(b)
2398                - d_i.instr_irrat(b) - d_i.incoherence(b);
2399       };
2400 
2401     simpler := cachedProcedure() {
2402       return procedure(d1, d2) {
2403                return k(d1.to_s) < k(d2.to_s);
2404              };
2405     };
2406 
2407     new := procedure() {
2408       return decision_algorithm({}, {}, {}, {}, {}, {}, {}, [], {}, {}, {},
2409                                 om, om, om);
2410     };
2411 
2412   }
2413 }
2414 
2415 // alias 
2416 da := procedure() { return decision_algorithm; };
   1 print("\nTesting Decision Algorithm");
   2 print(  "==========================");
   3 
   4 test_s := procedure() {
   5   d := decision_algorithm({'I'}, {'O'}, { ['A', 'P(A)'] }, // i, o, p
   6     { [['A', 'A'], 'P(A|A)'] }, {['A','U(A)']}, {'M'}, {'M2'}, // cp, u, m, m2
   7     [{ ['e', {['B','E_n1(B)']}], ['m2', {'M2_n1'}], ['o', {'O_n1'}] }], // n
   8     {}, {}, {}, om, om, om); // r, t_s, t_o, cacheds
   9   assert(d.s() == { 'P(A)', 'P(A|A)', 'U(A)', 'M', 'M2',
  10                     'E_n1(B)', 'M2_n1', 'O_n1' },
  11          "Error: s is incorrect.");
  12   print('[ok] s');
  13 };
  14 test_s();
  15 
  16 test_t := procedure() {
  17   d := decision_algorithm.new();
  18   d.i := { 'I' };
  19   d.o := { 'O' };
  20   d.p := { ['A', 'P'] };
  21   d.t_s := { // is                         -> s
  22              [{ ['I', false], ['P', false] }, { ['P', false] }],
  23              [{ ['I', false], ['P', true ] }, { ['P', true ] }],
  24              [{ ['I', true ], ['P', false] }, { ['P', true ] }],
  25              [{ ['I', true ], ['P', true ] }, { ['P', false] }]
  26            }; // xor
  27   d.t_o := { // p u m2       -> o
  28              [{ ['P', false] }, { ['O', false] }],
  29              [{ ['P', true ] }, { ['O', true ] }]
  30            }; // identity
  31   assert(d.t() == {
  32            // is                         ->  so
  33            [{["I", false], ["P", false]}, {["O", false], ["P", false]}],
  34            [{["I", false], ["P", true ]}, {["O", true ], ["P", true ]}],
  35            [{["I", true ], ["P", false]}, {["O", false], ["P", true ]}],
  36            [{["I", true ], ["P", true ]}, {["O", true ], ["P", false]}]
  37          },
  38          "Error: t is incorrect.");
  39   print('[ok] t');
  40 };
  41 test_t();
  42 
  43 test_ns := procedure() {
  44   d := decision_algorithm.new();
  45   d.n := [{ ['e', {['A','E_n1(A)']}], ['m2', {'M2_n1'}], ['o', {'O_n1'}] },
  46           { ['e', {['B','E_n2(B)']}], ['m2', {'M2_n2'}], ['o', {'O_n2'}] } ];
  47   assert(d.ns( 'e') == { 'E_n1(A)', 'E_n2(B)' } &&
  48          d.ns('m2') == {   'M2_n1', 'M2_n2'   } &&
  49          d.ns( 'o') == {    'O_n1', 'O_n2'    },
  50          "Error: ns is incorrect.");
  51   print('[ok] ns');
  52 };
  53 test_ns();
  54 
  55 test_new := procedure() {
  56   assert(decision_algorithm.new().t_o == {},
  57          "Error: new is incorrect.");
  58   print('[ok] new');
  59 };
  60 test_new();
  61 
  62 test_p_u_vars := procedure() {
  63   d := decision_algorithm.new();
  64   d.p  := { [       'A',   'P(A)'] };
  65   d.cp := { [['A', 'A'], 'P(A|A)'] };
  66   d.u  := { [       'A',   'U(A)'] };
  67   assert(d.p_u_vars() == {'P(A)', 'P(A|A)', 'U(A)'},
  68          "Error: is incorrect.");
  69   print('[ok] p_u_vars');
  70 };
  71 test_p_u_vars();
  72 test_p_vars := procedure() { test_p_u_vars(); };
  73 
  74 test_just_p_u_m2 := procedure() {
  75   d := decision_algorithm.new();
  76   d.i  := { 'I' };
  77   d.p  := { ['A', 'P(A)'] };
  78   d.u  := { ['A', 'U(A)'] };
  79   d.m2 := { 'M2' };
  80   assert(d.just_p_u_m2({ ['I', true],
  81                          ['P(A)', 0.3], ['U(A)', 4], ['M2', true] }) ==
  82                        { ['P(A)', 0.3], ['U(A)', 4], ['M2', true] },
  83          "Error: just_p_u_m2 is incorrect.");
  84   print('[ok] just_p_u_m2');
  85 };
  86 test_just_p_u_m2();
  87 
  88 test_parents := procedure() {
  89   d := decision_algorithm.new();
  90   d.m := { 'M1', 'M2' };
  91   d.r := { ['M1', {true, false}], ['M2', {true,false}] };
  92   d.t_s := { // is                  -> s
  93     [{ ['M1', false], ['M2', false] }, { ['M1', true ], ['M2', true ] }],
  94     [{ ['M1', false], ['M2', true ] }, { ['M1', false], ['M2', true ] }],
  95     [{ ['M1', true ], ['M2', false] }, { ['M1', false], ['M2', false] }],
  96     [{ ['M1', true ], ['M2', true ] }, { ['M1', true ], ['M2', false] }]
  97   }; // M1 = ~xor; M2 = ~M1
  98   assert(d.parents('M1') == {'M1', 'M2'} && d.parents('M2') == {'M1'},
  99          "Error: parents is incorrect.");
 100   print('[ok] parents');
 101   assert(d.children('M1') == { 'M1', 'M2' } && d.children('M2') == {'M1'},
 102         "Error: children is incorrect.");
 103   print('[ok] children');
 104 };
 105 test_children               := procedure() { test_parents(); };
 106 test_is_necessary_parent_of := procedure() { test_parents(); };
 107 test_vals_diff              := procedure() { test_parents(); };
 108 test_parents();
 109 
 110 /*   ,------.
 111     I--P--Y--X
 112      `--U-'
 113          `--Z
 114 */
 115 test_non_shared_ancestors := procedure() {
 116   d := decision_algorithm.new();
 117   d.i := {'I'};
 118   d.p := {['A', 'P']};
 119   d.u := {['B', 'U']};
 120   d.m2 := {'X', 'Y', 'Z'};
 121   d.r := { ['I', tf()], ['P', {0,1}], ['U', {0,1}],
 122            ['X', tf()], ['Y', tf()], ['Z', tf()] };
 123   ts_proc := procedure(is_state) {
 124     return { ['P', !is_state['I']], ['U', is_state['I']],
 125              ['X', !is_state['Y'] || is_state['I']],
 126              ['Y', is_state['P'] + is_state['U'] % 2],
 127              ['Z', is_state['U'] == 0] };
 128   };
 129   d.t_s := proc2func(ts_proc, d.poss_is_states());
 130   assert(d.non_shared_ancestors('X', d.m2) == { 'I', 'P', 'U' } &&
 131          d.non_shared_ancestors('Y', d.m2) == { 'P', 'U' } &&
 132          d.non_shared_ancestors('Z', d.m2) == { 'U' },
 133          "Error: non_shared_ancestors is incorrect.");
 134   print('[ok] non_shared_ancestors');
 135 
 136   assert(d.poss_i_states() == { {['I', true]}, {['I', false]} },
 137          "Error: poss_i_states is incorrect.");
 138   print('[ok] poss_i_states');
 139 
 140   assert({ {["P", 0], ["U", 0], ["X",  true], ["Y",  true], ["Z",  true]},
 141            {["P", 1], ["U", 1], ["X", false], ["Y", false], ["Z", false]}
 142          } < d.poss_s_states(),
 143          "Error: poss_s_states is incorrect");
 144   print('[ok] poss_s_states');
 145 
 146   assert(d.poss_u_states() == { {['U', 0]}, {['U', 1]} },
 147          "Error: poss_u_states is incorrect.");
 148   print('[ok] poss_u_states');
 149 
 150   assert(d.is_equal_to(d),
 151          "Error: is_equal_to is incorrect");
 152   print('[ok] is_equal_to true');
 153 
 154   d2 := d;
 155   d2.i := d.i + { 'I2' };
 156   assert(!d2.is_equal_to(d),
 157          "Error: is_equal_to is incorrect");
 158   print('[ok] is_equal_to false');
 159 };
 160 test_non_shared_ancestors();
 161 test_poss_is_states := procedure() { test_non_shared_ancestors(); };
 162 test_poss_i_states  := procedure() { test_non_shared_ancestors(); };
 163 test_poss_s_states  := procedure() { test_non_shared_ancestors(); };
 164 test_poss_u_states  := procedure() { test_non_shared_ancestors(); };
 165 test_is_equal_to    := procedure() { test_non_shared_ancestors(); };
 166 
 167 
 168 
 169 /* X--Y--O
 170     `--M--P
 171 */
 172 test_non_shared_descendants := procedure() {
 173   d := decision_algorithm.new();
 174   d.m2 := {'X', 'Y'};
 175   d.o := {'O'};
 176   d.m := {'M'};
 177   d.p := {['A','P']};
 178   d.r := { ['X', tf()], ['Y', tf()], ['O', tf()], ['M', tf()], ['P', {0,1}] };
 179   ts_proc := procedure(is_state) {
 180     return { ['X', true], ['Y', !is_state['X']], ['O', is_state['Y']],
 181              ['M', is_state['X']], ['P', #{ j : j in {1} | is_state['M']}] };
 182   };
 183   d.t_s := proc2func(ts_proc, d.poss_is_states());
 184   to_proc := procedure(p_u_m2) {
 185     return {[ 'O', p_u_m2['Y'] ]};
 186   };
 187   d.t_o := proc2func(to_proc, { d.just_p_u_m2(is) : is in d.poss_is_states() });
 188   assert(d.non_shared_descendants('X', d.m2) == {'M', 'O'} &&
 189          d.non_shared_descendants('Y', d.m2) == {'O'},
 190          "Error: non_shared_descendants is incorrect.");
 191   print('[ok] non_shared_descendants');
 192 };
 193 test_non_shared_descendants();
 194 
 195 d1 := cachedProcedure() {
 196   d   := decision_algorithm.new();
 197   d.i := { 'I' };
 198   d.p := { ['A', 'P'] };
 199   d.o := { 'O' };
 200   d.r := { [var, tf()] : var in {'I','O'} } + { ['P', {0,1}] };
 201   d.t_s := { [{['I', false], ['P', 0]}, {['P', 1]}], // eg P_t2 = ~I xor P_t1
 202              [{['I', false], ['P', 1]}, {['P', 0]}],
 203              [{['I',  true], ['P', 0]}, {['P', 0]}],
 204              [{['I',  true], ['P', 1]}, {['P', 1]}]
 205            };
 206   d.t_o := { [ {['P', 0]}, {['O',  true]} ],  // O = ~P
 207              [ {['P', 1]}, {['O', false]} ] };
 208   return d;
 209 };
 210 
 211 b1 := cachedProcedure() {
 212   b   := causal_markov_model.new();
 213   b.u := { 'U' };
 214   b.v := { 'V_P1', 'V_P2', 'O' };
 215   b.r := { [var, tf()] : var in {'V_P1', 'V_P2', 'O'} } + {['U',{1,2,3}]};
 216   /* P = V_P1 <-> V_P2
 217      O = ~ P
 218      O = ~ (V_P1 <-> V_P2) 
 219   */
 220   b.f['O'] := proc2func(procedure(cmm_st) {
 221                           return !(cmm_st['V_P1'] == cmm_st['V_P2']);
 222                         }, b.poss_states());
 223   /* P_t2 = I xor P_t1
 224           = U xor (V_P1 <-> V_P2)     
 225     t1  U V_P1 V_P2 P  ~I  t2  P  V_P1  V_P2
 226         1   F   F   T   T      F   T    F
 227         1   F   T   F   T      T   T    T
 228         1   T   F   F   T      T   F    F
 229         1   T   T   T   T      F   F    T
 230         2   F   F   T   T      F    
 231         2   F   T   F   T      T
 232         2   T   F   F   T      T
 233         2   T   T   T   T      F
 234         3   F   F   T   F      T   F    F    
 235         3   F   T   F   F      F   T    F
 236         3   T   F   F   F      F   F    T
 237         3   T   T   T   F      T   T    T
 238 
 239     two ways for U to be false {1,2}. {3} is truth
 240   */
 241   b.f['V_P1'] := proc2func(procedure(cmm_st) {
 242                    if (cmm_st['U'] == 3) {
 243                      return cmm_st['V_P2'];
 244                    } else {
 245                      return !cmm_st['V_P1'];
 246                    }
 247                  }, b.poss_states());
 248   b.f['V_P2'] := proc2func(procedure(cmm_st) {
 249                    if (cmm_st['U'] == 3) {
 250                      return cmm_st['V_P1'];
 251                    } else {
 252                      return cmm_st['V_P2'];
 253                    }
 254                  }, b.poss_states());
 255   b.n := 3;
 256   return b;
 257 };
 258 
 259 plf1 := cachedProcedure() {
 260   return { [{ ['U', 1]                         }, { ['I', false] }],
 261            [{ ['U', 2]                         }, { ['I', false] }],
 262            [{ ['U', 3]                         }, { ['I',  true] }],
 263            [{ ['V_P1', false], ['V_P2', false] }, { ['P',     1] }],
 264            [{ ['V_P1', false], ['V_P2',  true] }, { ['P',     0] }],
 265            [{ ['V_P1',  true], ['V_P2', false] }, { ['P',     0] }],
 266            [{ ['V_P1',  true], ['V_P2',  true] }, { ['P',     1] }],
 267            [{ ['O', false]                     }, { ['O', false] }],
 268            [{ ['O',  true]                     }, { ['O',  true] }]
 269          };
 270 };
 271 
 272 pf1 := cachedProcedure() {
 273   return d1().plf_to_pf(b1(), plf1());
 274 };
 275 
 276 test_poss_states := procedure() {
 277   assert(d1().poss_states() == {
 278            {["I", false], ["O", false], ["P", 0]},
 279            {["I", false], ["O", false], ["P", 1]},
 280            {["I", false], ["O",  true], ["P", 0]},
 281            {["I", false], ["O",  true], ["P", 1]},
 282            {["I",  true], ["O", false], ["P", 0]},
 283            {["I",  true], ["O", false], ["P", 1]},
 284            {["I",  true], ["O",  true], ["P", 0]},
 285            {["I",  true], ["O",  true], ["P", 1]}
 286          }, "Error: poss_states is incorrect.");
 287   print('[ok] poss_states');
 288 };
 289 test_poss_states();
 290 
 291 test_commutes := procedure() {
 292   // print('pf ' + pf1());  
 293 
 294   // This may take a while so it's commented out but it worked last time.
 295   // assert(d1().commutes(pf1(), b1()) == true, "Error: commutes true is incorrect.");
 296   // print('[ok] commutes true');
 297   print('[..] commutes true');
 298 };
 299 test_commutes();
 300 
 301 test_lf_commutes := procedure() {
 302   d := d1();
 303   d.cached_f := pf1();
 304   assert(d.lf_commutes(plf1(), b1()) == true,
 305          "Error: lf_commutes true is incorrect");
 306   print('[ok] lf_commutes true');
 307 
 308   plf := plf1();
 309   plf[ {['O', false]} ] := { ['O', true] };
 310   assert(d.lf_commutes(plf, b1()) == false,
 311          "Error: lf_commutes false is incorrect");
 312   print('[ok] lf_commutes false');
 313 };
 314 test_lf_commutes();
 315 
 316 d2 := procedure() {
 317   d := d1();
 318   d.cached_f  := pf1();
 319   d.cached_lf := plf1();
 320   return d;
 321 };
 322 
 323 test_is_implemented_by := procedure() {
 324   assert(d2().is_implemented_by(b1()) == true,
 325          "Error: is_implemented_by is incorrect.");
 326   print('[ok] is_implemented_by');
 327 };
 328 test_is_implemented_by();
 329 
 330 test_lf_inv := procedure() {
 331   assert(d2().lf_inv(b1())[ {['P',1]} ] == {
 332            { ['V_P1', false], ['V_P2', false] },
 333            { ['V_P1',  true], ['V_P2',  true] }
 334          }, "Error: lf_inv is incorrect");
 335   print('[ok] lf_inv');
 336 };
 337 test_lf_inv();
 338 test_plf_to_pf := procedure() { test_lf_inv(); };
 339 
 340 test_lf_inv_at := procedure() {
 341   assert(d2().lf_inv_at(b1(), 2)[ {['P',1]} ] == {
 342            { ['V_P1_t2', false], ['V_P2_t2', false] },
 343            { ['V_P1_t2',  true], ['V_P2_t2',  true] }
 344          }, "Error: lf_inv_at is incorrect");
 345   print('[ok] lf_inv_at');
 346 };
 347 test_lf_inv_at();
 348 
 349 test_lf_cm := procedure() {
 350   assert(d2().lf_cm(b1(), {['V_P1_t2',false],['V_P2_t2',false]}) == {['P',1]},
 351          "Error: lf_cm is incorrect");
 352   print('[ok] lf_cm');
 353 };
 354 test_lf_cm();
 355 
 356 test_lf_cm_vs := procedure() {
 357   assert(d2().lf_cm_vs(b1(), {['V_P1_t2',false],['V_P2_t2',false] }) == {'P'},
 358          "Error: lf_cm_vs is incorrect");
 359   print('[ok] lf_cm_vs');
 360 };
 361 test_lf_cm_vs();
 362 
 363 test_lf_cm_p_inv := procedure() {
 364   assert(d2().lf_cm_p_inv(b, { ['V_P1_t2', false], ['V_P2_t2', false] }) == 'A',
 365          "Error: lf_cm_p_inv is incorrect");
 366   print('[ok] p_inv');
 367   print('[ok] lf_cm_p_inv');
 368 };
 369 test_lf_cm_p_inv();
 370 test_p_inv := procedure() { test_lf_cm_p_inv(); };
 371 
 372 test_commutes_f := procedure() {
 373   // print('b.f[O] ' + b.f['O']);
 374   b := b1();
 375   b.f['O'] := proc2func(procedure(cmm_st) {
 376                           return cmm_st['V_P1'] && cmm_st['V_P2'];
 377                         }, b.poss_states());
 378   // print('b.f[O] ' + b.f['O']);
 379   // bs := {["O", true], ["U", 3], ["V_P1", true], ["V_P2", true]};
 380   // print('pf[bs] ' + pf[bs]);
 381   // print('t()[pf[bs]] ' + d.t()[pf[bs]]);
 382   // print('b.resp ' + b.response(b.v, bs));
 383   // print('pf[b.resp] ' + pf[b.response(b.v, bs)]);
 384 
 385   // assert(d1().commutes(pf1(), b) == false, "Error: commutes false is incorrect.");
 386   // print('[ok] commutes false');
 387   print('[..] commutes false');
 388 };
 389 test_commutes_f();
 390 
 391 test_drop_i := procedure() {
 392   d := decision_algorithm.new();
 393   d.i := { 'P(I_1)', 'P(I_2)' };
 394   ds  := { ['P(A)', 1], ['P(I_1)', true], ['P(I_2)', false], ['M1', true] };
 395   assert(d.drop_i(ds) == { ['P(A)', 1], ['M1', true] },
 396          "Error: drop_i is incorrect.");
 397   print('[ok] drop_i');
 398 };
 399 test_drop_i();
 400 
 401 /*   
 402     I--P--Y--X--O
 403      `--U-'
 404          `--Z--W
 405 */
 406 test_m2s_are_valid := procedure() {
 407   d    := decision_algorithm.new();
 408   d.i  := { 'I' };
 409   d.p  := { ['A', 'P'] };
 410   d.u  := { ['B', 'U'] };
 411   d.m2 := { 'X', 'Y', 'Z' };
 412   d.o  := { 'W', 'O' };
 413   d.r  := { [var, tf()] : var in d.i + d.s() };
 414   ts_proc := procedure(is_state) {
 415     return { ['P', !is_state['I']], ['U', is_state['I']],
 416              ['X', !is_state['Y']],
 417              ['Y', is_state['P'] && is_state['U']],
 418              ['Z', is_state['U'] == 0] };
 419   };
 420   d.t_s := proc2func(ts_proc, d.poss_is_states());
 421   to_proc := procedure(p_u_m2) {
 422     return { ['O', p_u_m2['X']], ['W', p_u_m2['W']] };
 423   };
 424   d.t_o := proc2func(to_proc, { d.just_p_u_m2(is) : is in d.poss_is_states() });
 425   assert(d.m2s_are_valid() == true, "Error: is incorrect.");
 426   print('[ok] m2s_are_valid true');
 427 
 428   // I------X
 429   ts_proc2 := procedure(is_state) {
 430     return { ['P', !is_state['I']], ['U', is_state['I']],
 431              ['X', !is_state['Y'] || is_state['I']],
 432              ['Y', is_state['P'] && is_state['U']],
 433              ['Z', is_state['U'] == 0] };
 434   };
 435   d.t_s := proc2func(ts_proc2, d.poss_is_states());
 436   assert(d.m2s_are_valid() == false, "Error: is incorrect.");
 437   print('[ok] m2s_are_valid false');
 438 
 439   // W-----U
 440   ts_proc3 := procedure(is_state) {
 441     return { ['P', !is_state['I']], ['U', is_state['I'] && is_state['W']],
 442              ['X', !is_state['Y']],
 443              ['Y', is_state['P'] && is_state['U']],
 444              ['Z', is_state['U'] == 0] };
 445   };
 446   d.t_s := proc2func(ts_proc2, d.poss_is_states());
 447   assert(d.m2s_are_valid() == false, "Error: is incorrect.");
 448   print('[ok] m2s_are_valid false 2');
 449 };
 450 test_m2s_are_valid();
 451 
 452 /*
 453 test_no_e_in_u := procedure() {
 454   d := decision_algorithm.new();
 455   assert(d.no_e_in_u() == ,
 456          "Error: no_e_in_u is incorrect.");
 457   print('[ok] no_e_in_u false');
 458 
 459   assert(d.no_e_in_u() == ,
 460          "Error: no_e_in_u is incorrect.");
 461   print('[ok] no_e_in_u true');
 462 };
 463 test_no_e_in_u();
 464 */
 465 
 466 /* test_i_ref := procedure() {
 467   print(d1().i_ref(b2(), b2())); // out of memory
 468 //  assert(d1().i_ref(b1(), b1()) == 
 469 //         "Error: i_ref is incorrect.");
 470   print('[  ] i_ref');
 471 };
 472 test_i_ref();
 473 */
 474 
 475 test_str_exprs := procedure() {
 476   d := decision_algorithm.new();
 477   d.i := { 'P(I_1)' };
 478   d.p := { ['A', 'P(A)'], ['B', 'P(B)'], ['I_1', 'P(I_1)'],
 479            [And('A', 'B'), 'P(A^B)'],
 480            [BoxArrow('A', 'B'), 'P(A[]->B)'] };
 481   assert(d.str_exprs() == { ['A',0] , ['B',0] },
 482          "Error: is incorrect.");
 483   print('[ok] str_exprs');
 484 };
 485 test_str_exprs();
 486 
 487 test_preds := procedure() {
 488   d := decision_algorithm.new();
 489   d.i := { 'P(I_1)' };
 490   d.p := { ['A', 'P(A)'], ['B', 'P(B)'], ['I_1', 'P(I_1)'],
 491            [And('A', 'B'), 'P(A^B)'], [Q('A'), 'P(Q(A))'],
 492            [R('A','B'), 'P(R(A,B))'],
 493            [BoxArrow('A', 'B'), 'P(A[]->B)'] };
 494   assert(d.preds() == { ['Q', 1], ['R', 2] },
 495          "Error: is incorrect.");
 496   print('[ok] preds arity');
 497 };
 498 test_preds();
 499 
 500 w1 := procedure() {
 501   w := causal_markov_model.new();
 502   w.u := {};
 503   w.v := { 'X', 'Y' };
 504   w.n := 2;
 505   w.r := { ['X', tf()], ['Y', tf()] };
 506   w.f := { ['X', { [{['X', false], ['Y', false]}, false],
 507                    [{['X', false], ['Y',  true]},  true],
 508                    [{['X',  true], ['Y', false]},  true],
 509                    [{['X',  true], ['Y',  true]},  true]
 510                  }],
 511            ['Y', { [{['X', false], ['Y', false]},  true],
 512                    [{['X', false], ['Y',  true]},  true],
 513                    [{['X',  true], ['Y', false]},  true],
 514                    [{['X',  true], ['Y',  true]}, false]
 515                  }]
 516          };
 517 };
 518 
 519 test_poss_base_expr_refs := procedure() {
 520   d := decision_algorithm.new();
 521   w := w1();
 522   // print(d.poss_base_expr_refs(w, { ['A', 0] }));
 523   // print(d.poss_base_expr_refs(w, { ['A', 0], [Q('A'), 1], [R('A','B'), 2] }));
 524   // print(d.poss_base_expr_refs(w, { ['A', 0], [Q('A'), 1] }));
 525   // 2 ** 81 subsets
 526   // assert(d.poss_base_expr_refs(w, { 'A', Q('A'), R('A','B') }),
 527   //        "Error: is incorrect.");
 528   // print('[ok] poss_base_expr_refs');
 529 };
 530 test_poss_base_expr_refs();
 531 
 532 /*
 533 test_poss_base_refs := procedure() {
 534   d := decision_algorithm.new();
 535   d.cached_lf := { };
 536   w := w1();
 537   w_sts := w1.cm().poss_states();
 538   pbers := { 
 539     { ['A', { w_sts[1], w_sts[2] }],
 540       ['Q', { [{ w_sts[1], w_sts[2] }], 
 541               [{ w_sts[4] }] }] },
 542     { ['A', { w_sts[3], w_sts[4] }],
 543       ['Q', { [{ w_sts[2], w_sts[3] }], 
 544               [{ w_sts[5] }] }] }
 545   };  
 546   // assert(d.poss_base_refs(w, b, pbers) == ,
 547   //         "Error: poss_base_refs is incorrect.");
 548   print('[ok] poss_base_refs');
 549 };
 550 // test_poss_base_refs();
 551 */
 552 
 553 test_rr := procedure() {
 554   d := decision_algorithm.new();
 555   b := causal_markov_model.new();
 556   b.v := { 'V_1' };
 557   b.r := { ['V_1', tf()] };
 558   b.f := {
 559     ['V_1', { [{ ['V_1', false] },  true],
 560               [{ ['V_1',  true] }, false]
 561             }]
 562   };
 563   b.n := 3; // bug when set to 2?
 564   // rr(expr, [w, b, atom_ref, ti, context])
 565   assert(d.rr('A', [b, b, {}, 1, { ['A', { {['V_1_t1', true]} }] }]) ==
 566                                          { {['V_1_t1', true]} },
 567          "Error: rr context is incorrect.");
 568   print('[ok] rr context');
 569 
 570   d.i := { 'P(I_1)' };
 571   d.p := { ['I_1','P(I_1)'] };
 572   b.a := [ {['V_1_t1',  true]},
 573            {['V_1_t2', false]},
 574            {['V_1_t3',  true]} ];
 575   d.cached_lf := {
 576     [ {['V_1', false]}, {['P(I_1)', false]} ],
 577     [ {['V_1',  true]}, {['P(I_1)',  true]} ]
 578   };
 579   all_sts := b.cm().poss_states();
 580   assert(d.rr('I_1', [b, b, {}, 1, {}]) ==
 581            compat_compl_sts( {['V_1_t1', true]}, all_sts ),
 582          "Error: rr i is incorrect.");
 583   print('[ok] rr i');
 584 
 585   assert(d.rr(At_t('I_1', -1), [b, b, {}, 2, {}]) ==
 586            compat_compl_sts( {['V_1_t1', true]}, all_sts ),
 587          "Error: rr At_t is incorrect.");
 588   print('[ok] rr At_t');
 589 
 590   atom_ref := { ['B', compat_compl_sts( {['V_1_t1', true]}, all_sts )],
 591                 ['IsAorB', { [compat_compl_sts({['V_1_t1', true]}, all_sts)],
 592                              [compat_compl_sts({['V_1_t2', true]}, all_sts)] }]
 593               };
 594   assert(d.rr('B', [b, b, atom_ref, 1, {}]) ==
 595            compat_compl_sts({['V_1_t1', true]}, all_sts) &&
 596          d.rr('IsAorB', [b, b, atom_ref, 1, {}]) == {
 597            [compat_compl_sts({['V_1_t1', true]}, all_sts)],
 598            [compat_compl_sts({['V_1_t2', true]}, all_sts)]
 599          },
 600          "Error: rr atom_ref is incorrect.");
 601   print('[ok] rr atom_ref');
 602 
 603   atom_ref := { ['C', compat_compl_sts( {['V_1_t1', true]}, all_sts )],
 604                 ['D', compat_compl_sts( {['V_1_t2', true]}, all_sts )] };
 605   assert(d.rr(And('C', 'D'), [b, b, atom_ref, 1, {}]) ==
 606            compat_compl_sts({['V_1_t1', true], ['V_1_t2', true]}, all_sts),
 607          "Error: rr And is incorrect");
 608   print('[ok] rr And');
 609 
 610   assert(d.rr(Or('C', 'D'), [b, b, atom_ref, 1, {}]) ==
 611            compat_compl_sts({['V_1_t1', true]}, all_sts) +
 612            compat_compl_sts({['V_1_t2', true]}, all_sts),
 613          "Error: rr Or is incorrect");
 614   print('[ok] rr Or');
 615 
 616   assert(d.rr(Not('C'), [b, b, atom_ref, 1, {}]) ==
 617            compat_compl_sts({['V_1_t1', false]}, all_sts),
 618          "Error: rr Not is incorrect");
 619   print('[ok] rr Not');
 620 
 621   assert(d.rr(Implies('C', 'D'), [b, b, atom_ref, 1, {}]) ==
 622            compat_compl_sts({['V_1_t1', false]}, all_sts) +
 623            compat_compl_sts({['V_1_t2', true]}, all_sts),
 624          "Error: rr Implies is incorrect");
 625   print('[ok] rr Implies');
 626 
 627   print('[  ] rr');
 628 };
 629 test_rr();
 630 
 631 test_true_p := procedure() {
 632   d := decision_algorithm.new();
 633   b_event := {['V', true]};
 634   ref_i := { [ b_event, { {['X', true]}, {['Y', true]} }] };
 635   w := causal_markov_model.new();
 636   // TODO: rewrite to use states rather than events
 637 // or just fold into sq_err
 638   w.a := [{ ['Y', true] }];
 639   assert(d.true_p(w, ref_i, b_event) == 1, "Error: true_p true is incorrect.");
 640   print('[ok] true_p true');
 641 
 642   w.a := [{ ['X', true] }];
 643   assert(d.true_p(w, ref_i, b_event) == 1, "Error: true_p true is incorrect.");
 644   print('[ok] true_p true 2');
 645 
 646   w.a := [{ ['X', false] }];
 647   assert(d.true_p(w, ref_i, b_event) == 0, "Error: true_p false is incorrect.");
 648   print('[ok] true_p false');
 649 };
 650 //test_true_p();
 651 
 652 test_sq_err := procedure() {
 653   d := decision_algorithm.new();
 654   d.p := { ['A', 'P(A)'] };
 655   d.cached_lf := { [ {['V',  true]}, {['P(A)', 0.75]} ],
 656                    [ {['V', false]}, {['P(A)', 0.33]} ] };
 657   b := causal_markov_model.new();
 658   b.a := [{['V_t1', true]}, {['V_t2', false]}];
 659   ref_i := { [{['V_t1',  true]}, { {['V_t1', true], ['V_t2', false]},
 660                                    {['V_t1', true], ['V_t2',  true]} }],
 661              [{['V_t2', false]}, { {['V_t1', true], ['V_t2', false]},
 662                                    {['V_t1', true], ['V_t2',  true]} }] };
 663   assert(abs(d.sq_err(b, b, ref_i) - 0.5114) < 0.01,
 664          "Error: sq_err is incorrect.");
 665   print('[ok] true_p');
 666   print('[ok] sq_err');
 667 };
 668 test_sq_err();
 669 test_true_p := procedure() { test_sq_err(); };
 670 
 671 b2 := procedure() {
 672   b := b1();
 673   b.n := 2;
 674   b.a := [
 675     { ['U_t1', 1], ['V_P1_t1', false], ['V_P2_t1', false], ['O_t1',  true] },
 676     { ['U_t2', 2], ['V_P1_t2',  true], ['V_P2_t2', false], ['O_t2', false] }
 677   ];
 678   return b;
 679 };
 680 
 681 test_ref_expr := procedure() {
 682   d := d2();
 683   b := b2();
 684   bs1 := {{['U', 1], ['V_P1_t1', false], ['V_P2_t1', false], ['O', false]}};
 685   d.cached_ref := {
 686     [ {['V_P1_t1', false], ['V_P2_t1', false]}, bs1 ]
 687   };
 688   assert(d.ref_expr(b, b, 1, 'A') == bs1,
 689          "Error: ref_expr is incorrect.");
 690   print('[ok] ref_expr');
 691 };
 692 test_ref_expr();
 693 
 694 test_actual_p_events := procedure() {
 695   assert(d2().actual_p_events(b2()) == {
 696            { ['V_P1_t1', false], ['V_P2_t1', false] },
 697            { ['V_P1_t2',  true], ['V_P2_t2', false] }
 698          },
 699          "Error: actual_p_events is incorrect.");
 700   print('[ok] actual_p_events');
 701 };
 702 test_actual_p_events();
 703 
 704 test_prob_states := procedure() {
 705   assert(d2().prob_states(b2()) == { [1, {['P', 1]}], [2, {['P', 0]}] },
 706          "Error: prob_states is incorrect.");
 707   print('[ok] prob_states');
 708   print('[ok] prob_states_t');
 709 };
 710 test_prob_states();
 711 test_prob_states_t := procedure() { test_prob_states(); };
 712 
 713 d3 := procedure() {
 714   d := decision_algorithm.new();
 715   d.u := { ['A', 'U(A)'] };
 716   d.o := { 'O' };
 717   d.n := [{ ['e', {['B','E1(B)']}], ['m2',{'M2_1'}], ['o',{'O_1'}] },
 718           { ['e', {['C','E2(C)']}], ['m2',{'M2_2'}], ['o',{'O_2'}] }];
 719   return d;
 720 };
 721 
 722 test_eo := procedure() {
 723   assert(d3().eo() == { [{ ['A',  'U(A)'] }, { 'O' }],
 724                         [{ ['B', 'E1(B)'] }, {'O_1'}],
 725                         [{ ['C', 'E2(C)'] }, {'O_2'}] },
 726          "Error: eo is incorrect.");
 727   print('[ok] eo');
 728 };
 729 test_eo();
 730 
 731 test_instr_irrat := procedure() {
 732   print('test instr irrat');
 733   d := decision_algorithm.new();
 734   d.o := { 'O' };
 735   d.p := { ['A', 'P(A)'],
 736            [Not('A'), 'P(~A)'],
 737            [BoxArrow('A', 'B'), 'P(A []-> B)'],
 738            [BoxArrow(Not('A'), 'B'), 'P(~A []-> B)'] };
 739   d.u := { ['B', 'U(B)'] };
 740   d.r := { ['P(A)', {0.3, 0.7}],
 741            ['P(~A)', {0.1, 0.9}],
 742            ['P(A []-> B)', {0.2, 0.8}],
 743            ['P(~A []-> B)', {0.25, 0.75}],
 744            ['U(B)', { 5 }], ['O', tf()] };
 745   b := causal_markov_model.new();
 746   b.u := {};
 747   b.v := { 'V_A', 'V_~A', 'V_A_[]->_B', 'V_~A_[]->_B', 'V_O', 'V_U' };
 748   b.r := { [var, tf()] : var in b.v };
 749   false_func := proc2func(procedure(bs) { return false; }, b.poss_states());
 750   print('false func');
 751   b.f := {
 752     ['V_A', false_func],
 753     ['V_~A', false_func],
 754     ['V_A_[]->_B', false_func],
 755     ['V_~A_[]->_B', false_func],
 756     ['V_O', false_func],
 757     ['V_U', false_func]
 758   };
 759   b.n := 3;
 760   b.a := [
 761     { ['V_A_t1',  true], ['V_~A_t1', false],
 762       ['V_~A_[]->_B_t1', false], ['V_~A_[]->_B_t1', true],
 763       ['V_O_t1', false], ['V_U_t1', false] },
 764     { ['V_A_t2',  true], ['V_~A_t2', false],
 765       ['V_~A_[]->_B_t2', false], ['V_~A_[]->_B_t2', true],
 766       ['V_O_t2', false], ['V_U_t2', false] },
 767     { ['V_A_t3',  true], ['V_~A_t3', false],
 768       ['V_~A_[]->_B_t3', false], ['V_~A_[]->_B_t3', true],
 769       ['V_O_t3', false], ['V_U_t3', false] }
 770   ];
 771   d.cached_lf := {
 772     [{['V_A', false]},         {['P(A)', 0.3]}        ],
 773     [{['V_A',  true]},         {['P(A)', 0.7]}        ],
 774     [{['V_~A', false]},        {['P(~A)', 0.1]}       ],
 775     [{['V_~A',  true]},        {['P(~A)', 0.9]}       ],
 776     [{['V_A_[]->_B', false]},  {['P(A []-> B)', 0.2]} ],
 777     [{['V_A_[]->_B',  true]},  {['P(A []-> B)', 0.8]} ],
 778     [{['V_~A_[]->_B', false]}, {['P(~A []-> B)', 0.2]}],
 779     [{['V_~A_[]->_B',  true]}, {['P(~A []-> B)', 0.8]}],
 780     [{['V_O', false]},         {['O', false]}         ],
 781     [{['V_O',  true]},         {['O',  true]}         ],
 782     [{['V_U', false]},         {['U(B)', false]}      ],
 783     [{['V_U',  true]},         {['U(B)',  true]}      ]
 784   };
 785   print('cached lf');
 786   d.cached_f := d.plf_to_pf(b, d.cached_lf);
 787   print('cached f');
 788   b_cm := b.cm();
 789   print('b cm');
 790   b_cm_poss_states := b_cm.poss_states();
 791   print(#b_cm_poss_states);
 792   o_t2 := compat_compl_sts({['V_O_t2', true]}, b_cm_poss_states);
 793   print('o_t2');
 794   o_t3 := compat_compl_sts({['V_O_t3', true]}, b_cm_poss_states);
 795   no_t2 := compat_compl_sts({['V_O_t2', false]}, b_cm_poss_states);
 796   no_t3 := compat_compl_sts({['V_O_t3', false]}, b_cm_poss_states);
 797   print('no_t3');
 798   d.cached_ref := {
 799     [{['V_A_t1', true ]}, o_t2],
 800     [{['V_A_t1', false]}, o_t2],
 801 
 802     [{['V_~A_t1', true ]}, no_t2],
 803     [{['V_~A_t1', false]}, no_t2],
 804 
 805     [{['V_A_t2', true ]}, o_t3],
 806     [{['V_A_t2', false]}, o_t3],
 807 
 808     [{['V_~A_t2', true ]}, no_t3],
 809     [{['V_~A_t2', false]}, no_t3]
 810 
 811   };
 812   // print(d.instr_irrat(b, b));
 813   // assert(d.instr_irrat() == ,
 814   //       "Error: instr_irrat is incorrect.");
 815   print('[.?] instr_irrat');
 816 };
 817 // test_instr_irrat();
 818 
 819 test_max_p_dist := procedure() {
 820   d := decision_algorithm.new();
 821   d.p := { ['A', 'P(A)'], ['B', 'P(B)'] };
 822   thirds := { 0, 0.33, 0.66, 1 };
 823   d.r := { ['P(A)', thirds], ['P(B)', thirds] };
 824   b := causal_markov_model.new();
 825   b.v := { 'V_A', 'V_B' };
 826   four := { 1,2,3,4 };
 827   b.r := { ['V_A', four], ['V_B', four] };
 828   d.cached_f := {
 829     [{['V_A', 1], ['V_B', 1]}, {['P(A)',    0], ['P(B)',    0]}],
 830     [{['V_A', 1], ['V_B', 2]}, {['P(A)',    0], ['P(B)', 0.33]}],
 831     [{['V_A', 1], ['V_B', 3]}, {['P(A)',    0], ['P(B)', 0.66]}],
 832     [{['V_A', 1], ['V_B', 4]}, {['P(A)',    0], ['P(B)',    1]}],
 833     [{['V_A', 2], ['V_B', 1]}, {['P(A)', 0.33], ['P(B)',    0]}],
 834     [{['V_A', 2], ['V_B', 2]}, {['P(A)', 0.33], ['P(B)', 0.33]}],
 835     [{['V_A', 2], ['V_B', 3]}, {['P(A)', 0.33], ['P(B)', 0.66]}],
 836     [{['V_A', 2], ['V_B', 4]}, {['P(A)', 0.33], ['P(B)',    1]}],
 837     [{['V_A', 3], ['V_B', 1]}, {['P(A)', 0.66], ['P(B)',    0]}],
 838     [{['V_A', 3], ['V_B', 2]}, {['P(A)', 0.66], ['P(B)', 0.33]}],
 839     [{['V_A', 3], ['V_B', 3]}, {['P(A)', 0.66], ['P(B)', 0.66]}],
 840     [{['V_A', 3], ['V_B', 4]}, {['P(A)', 0.66], ['P(B)',    1]}],
 841     [{['V_A', 4], ['V_B', 1]}, {['P(A)',    1], ['P(B)',    0]}],
 842     [{['V_A', 4], ['V_B', 2]}, {['P(A)',    1], ['P(B)', 0.33]}],
 843     [{['V_A', 4], ['V_B', 3]}, {['P(A)',    1], ['P(B)', 0.66]}],
 844     [{['V_A', 4], ['V_B', 4]}, {['P(A)',    1], ['P(B)',    1]}]
 845   };
 846   b.a := [{['V_A_t1', 2], ['V_B_t1', 3]}];
 847   assert(d.max_p_dist(b, 1) == 1.33,
 848          "Error: max_p_dist is incorrect.");
 849   print('[ok] max_p_dist');
 850 };
 851 test_max_p_dist();
 852 
 853 test_poss_cms := procedure() {
 854   d := decision_algorithm.new();
 855   d.p := { ['A', 'P(A)'], ['B', 'P(B)'], ['C', 'P(C)'] };
 856   d.r := { [var, tf()] : var in range(d.p) };
 857   assert(exists(cm in d.poss_cms() |
 858            cm.f == {
 859              ["P(A)", { [{["P(B)", false]},  true],
 860                         [{["P(B)",  true]}, false] }],
 861              ["P(B)", { [{["P(C)", false]}, false],
 862                         [{["P(C)",  true]},  true] }]
 863            }) &&
 864          exists(cm in d.poss_cms() |
 865            cm.f == {
 866              ["P(A)", { [{["P(B)", false]}, false],
 867                         [{["P(B)",  true]}, false]} ],
 868              ["P(C)", { [{["P(A)", false]},  true],
 869                         [{["P(A)",  true]},  true] }]
 870            }),
 871          "Error: poss_cms is incorrect.");
 872   print('[ok] poss_cms');
 873 };
 874 test_poss_cms();
 875 test_poss_cms_from_l      := procedure() { test_poss_cms(); };
 876 test_poss_cms_from_p_pars := procedure() { test_poss_cms(); };
 877 
 878 test_e_func := procedure() {
 879   d := d3();
 880   ds := { ['U(A)', 5], ['E1(B)', 3], ['E2(C)', 8] };
 881   assert(d.e_func(ds, 0) == { ['U(A)',  5] } &&
 882          d.e_func(ds, 1) == { ['E1(B)', 3] } &&
 883          d.e_func(ds, 2) == { ['E2(C)', 8] },
 884          "Error: e_func is incorrect.");
 885   print('[ok] e_func');
 886 };
 887 test_e_func();
 888 
 889 test_ue := procedure() {
 890   d := d3();
 891   assert(d.ue(0) == d.u && d.ue(1) == d.n[1]['e'] && d.ue(2) == d.n[2]['e'],
 892          "Error: ue is incorrect.");
 893   print('[ok] ue');
 894 };
 895 test_ue();
 896 
 897 test_ue_base_exprs := procedure() {
 898   d := d3();
 899   assert(d.ue_base_exprs(0) == { 'A' } &&
 900          d.ue_base_exprs(1) == { 'B' } &&
 901          d.ue_base_exprs(2) == { 'C' },
 902          "Error: ue_base_exprs is incorrect.");
 903   print('[ok] ue_base_exprs');
 904 };
 905 test_ue_base_exprs();
 906 
 907 /*
 908 test_state_space := procedure() {
 909   d := decision_algorithm.new();
 910   d.p := { ['A', 'P(A)'], ['B', 'P(B)'] };
 911   assert(d.state_space(range(d.p)) == { 
 912            {Not('A'), Not('B')},
 913            {Not('A'),     'B' },
 914            {    'A' , Not('B')},
 915            {    'A' ,     'B' }
 916          },
 917          "Error: state_space is incorrect.");
 918   print('[ok] state_space');
 919   assert(d.event_space(range(d.p)) == d.state_space(range(d.p)) + {
 920            { Not('A') }, { 'A' }, { Not('B') }, { 'B' }, {}
 921          },
 922          "Error: event_space is incorrect.");
 923   print('[ok] event_space');
 924   // print(d.expr_space(d.state_space(range(d.p))));
 925   assert(d.expr_space( d.state_space(range(d.p)) ) == {
 926            And("A",      "B"), And(    "A" , Not("B")), 
 927            And("B", Not("A")), And(Not("A"), Not("B"))
 928          },
 929          "Error: expr_space is incorrect.");
 930   print('[ok] expr_space');
 931 
 932   d.p[And("A", "B")] := 'P(A^B)';
 933   d.p[And("A", Not("B"))] := 'P(A^~B)';
 934   d.p[And("B", Not("A"))] := 'P(B^~A)';
 935   d.p[And(Not("A"), Not("B"))] := 'P(~A^~B)';
 936   d.r := { ['P(A^B)', {0,1}], ['P(A^~B)', {0,1}], 
 937            ['P(B^~A)', {0,1}], ['P(~A^~B)', {0,1}] };
 938   assert({["P(A^B)", 1], ["P(A^~B)", 0], ["P(B^~A)", 0], ["P(~A^~B)", 0]}
 939            in d.poss_p_vs({'P(A)', 'P(B)'}),
 940          "Error: poss_p_vs is incorrect");
 941   print('[ok] poss_p_vs');
 942 };
 943 test_state_space();
 944 test_event_space := procedure() { test_state_space(); };
 945 test_expr_space := procedure() { test_expr_space(); };
 946 test_poss_p_vs := procedure() { test_poss_p_vs(); };
 947 */
 948 
 949 // deprecated?
 950 /*
 951 test_p_of_vs := procedure() {
 952   d := decision_algorithm.new();
 953   d.p := { ['A', 'P(A)'], ['B', 'P(B)'], ['C', 'P(C)'] };
 954   cm := causal_model(u, v, r, f, om);
 955   cm.u := { 'P(A)', 'P(B)' };
 956   cm.v := { 'P(C)' };
 957   cm.r := { [var, tf()] : var in cm.u + cm.v };
 958   cm.f := { 
 959     ['P(C)', { [ {['P(A)', false], ['P(B)', false]}, false],
 960                [ {['P(A)', false], ['P(B)',  true]},  true],
 961                [ {['P(A)',  true], ['P(B)', false]},  true],
 962                [ {['P(A)',  true], ['P(B)',  true]},  true]
 963              }]
 964   };
 965   assert(d.p_of_vs(cm, { ['P(A)', 0.33], ['P(B)', 0.5] }) == ,
 966          "Error: p_of_vs is incorrect.");
 967   print('[ok] p_of_vs');
 968 };
 969 test_p_of_vs();
 970 */
 971 
 972 test_coh_pcms := procedure() {
 973   cm := causal_model.new();
 974   cm.u := { 'A' };
 975   cm.v := { 'B' };
 976   cm.r := { ['A', tf()], ['B', tf()] };
 977   cm.f := { ['B', { [{ ['A', false] },  true],
 978                     [{ ['A',  true] }, false] }] };
 979   d := decision_algorithm.new();
 980   d.p := { ['A', 'P(A)'] };
 981   d.r := { ['P(A)', { 0, 1 }] };
 982 //  print(d.coh_pcms(cm));
 983 //  assert(d.coh_pcms(cm) == ,
 984 //         "Error: coh_pcms is incorrect.");
 985   print('[.?] coh_pcms');
 986 };
 987 // test_coh_pcms();
 988 
 989 /* Template
 990 test_ := procedure() {
 991   d := decision_algorithm.new();
 992   assert(d.() == ,
 993          "Error:  is incorrect.");
 994   print('[ok] ');
 995 };
 996 test_();
 997 */
 998 
 999 print(' ');
1000 
1001 test_expr_to_p_funcs := procedure() {
1002   d := decision_algorithm.new();
1003   d.p := { ['A', 'P(A)'], ['B', 'P(B)'] };
1004   assert(d.expr_to_p_funcs(Implies('A', 'B'), {'A', 'B' }) == {
1005            {["P(A)", false], ["P(B)", false]},
1006            {["P(A)", false], ["P(B)",  true]},
1007            {["P(A)",  true], ["P(B)",  true]}
1008          }, "Error: expr_to_p_funcs is incorrect.");
1009   print('[ok] expr_to_p_funcs');
1010 
1011   d.p[Implies('A', 'B')] := 'P(A->B)';
1012   d.p[Not(And('A', Not('B')))] := 'P(~(A^~B))';
1013   assert(d.equiv_p_expr(d.p, d.cp, 'P(A->B)', 'P(~(A^~B))') == true,
1014          "Error: equiv_p_expr is incorrect");
1015   print('[ok] equiv_p_expr true');
1016 };
1017 test_expr_to_p_funcs();
1018 
1019 test_sum_p := procedure() {
1020   d := decision_algorithm.new();
1021   d.p := { ['A', 'P(A)'], ['B', 'P(B)'], ['C', 'P(C)'] };
1022   p1 := { ['P(A)', 0.6], ['P(C)', 0.3] };
1023   p2 := { ['P(B)', 0.2], ['P(C)', 0.5] };
1024   assert(d.sum_p(p1, p2) == { ['P(A)', 0.6], ['P(B)', 0.2], ['P(C)', 0.8] },
1025          "Error: sum_p is incorrect.");
1026   print('[ok] sum_p');
1027 };
1028 test_sum_p();
1029 
1030 test_add_p_distr := procedure() {
1031   d := decision_algorithm.new();
1032   d.p := { ['A', 'P(A)'], ['B', 'P(B)'], ['C', 'P(C)'] };
1033   d.cp := { [['C', 'A'], 'P(C|A)'], [['C', 'B'], 'P(C|B)'] };
1034   p_cp1 := [{ ['P(A)', 0.2] }, { ['P(C|A)', 0.3], ['P(C|B)', 0.4] }];
1035   p_cp2 := [{ ['P(A)', 0.4] }, { ['P(C|A)', 0.1], ['P(C|B)', 0.3] }];
1036   p_cp := d.add_p_distr(p_cp1, p_cp2);
1037   assert((p_cp[2])['P(C|A)'] == 0.4 && (p_cp[2])['P(C|B)'] == 0.7 &&
1038          (p_cp[1])['P(A)'] - 0.6 < 0.01,
1039          "Error: add_p_distr is incorrect.");
1040   print('[ok] add_p_distr');
1041 };
1042 test_add_p_distr();
1043 
1044 test_prob_distance := procedure() {
1045   d := decision_algorithm.new();
1046   d.p := { ['A', 'P(A)'], ['B', 'P(B)'] };
1047   d.cp := { [['B', 'A'], 'P(B|A)'] };
1048   pp1 := { ['P(A)', 0.4], ['P(B)', 0.5], ['P(B|A)', 0.6] };
1049   pp2 := { ['P(A)', 0.1], ['P(B)', 0.9], ['P(B|A)', 0.2] };
1050   assert(d.prob_distance(d.p, d.cp, pp1, pp2) == 1.1,
1051          "Error: prob_distance is incorrect.");
1052   print('[ok] prob_distance');
1053 };
1054 test_prob_distance();
1055 
1056 test_poss_io_states := procedure() {
1057   d := decision_algorithm.new();
1058   d.i := { 'P(I)' };
1059   d.o := { 'O' };
1060   d.r := { ['P(I)', {0,1}], ['O', tf()]};
1061   assert(d.poss_io_states() == {
1062           { ['P(I)', 0], ['O', false] },
1063           { ['P(I)', 0], ['O',  true] },
1064           { ['P(I)', 1], ['O', false] },
1065           { ['P(I)', 1], ['O',  true] }
1066          },
1067          "Error: poss_io_states is incorrect.");
1068   print('[ok] poss_io_states');
1069 };
1070 test_poss_io_states();
1071 
1072 test_poss_m_states := procedure() {
1073   d := decision_algorithm.new();
1074   d.m := { 'M1', 'M2' };
1075   d.r := { ['M1', tf()], ['M2', tf()] };
1076   assert(d.poss_m_states() == {
1077            { ['M1', false], ['M2', false] },
1078            { ['M1', false], ['M2',  true] },
1079            { ['M1',  true], ['M2', false] },
1080            { ['M1',  true], ['M2',  true] }
1081          },
1082          "Error: poss_m_states is incorrect.");
1083   print('[ok] poss_m_states');
1084 };
1085 test_poss_m_states();
1086 
1087 test_poss_o_states := procedure() {
1088   d := decision_algorithm.new();
1089   d.o := { 'O1', 'O2' };
1090   d.r := { ['O1', tf()], ['O2', tf()] };
1091   assert(d.poss_o_states() == {
1092            { ['O1', false], ['O2', false] },
1093            { ['O1', false], ['O2',  true] },
1094            { ['O1',  true], ['O2', false] },
1095            { ['O1',  true], ['O2',  true] }
1096          },
1097          "Error: poss_o_states is incorrect.");
1098   print('[ok] poss_o_states');
1099 };
1100 test_poss_o_states();
1101 
1102 test_possible_fms := procedure() {
1103   d := decision_algorithm.new();
1104   d.o := { 'O' };
1105   d.m := { 'M_1', 'M_2' };
1106   d.m2 := { 'M2_1', 'M2_2' };
1107   d.r := { [var, tf()] : var in d.o + d.m + d.m2 };
1108   // d.r := { [var, tf()] : var in d.o + d.m };
1109   ds := decision_algorithm.new();
1110   ds.o := { 'O' };
1111   ds.m := { 'ds.M_1', 'ds.M_2' };
1112   ds.m2 := { 'ds.M2_1', 'ds.M2_2' };
1113   ds.r := { [var, tf()] : var in ds.o + ds.m + ds.m2 };
1114   // ds.r := { [var, tf()] : var in ds.o + ds.m };
1115   p_fms := [fm : fm in d.possible_fms(ds)];
1116   print(p_fms);
1117   // print('[.?] possible_fms');
1118 };
1119 // test_possible_fms();
1120 test_apply_fm_m := procedure() { test_possible_fms(); };
1121 test_fm_commutes := procedure() { test_possible_fms(); };
1122 
1123 test_force_list := procedure() {
1124   d := decision_algorithm.new();
1125   d.p := { ['A', 'P(A)'] };
1126   d.cp := { [['B', 'A'], 'P(B|A)'] };
1127   assert(d.force_list('P(A)'  ) == ['A',  om] &&
1128          d.force_list('P(B|A)') == ['B', 'A'],
1129          "Error: force_list  is incorrect.");
1130   print('[ok] force_list ');
1131 };
1132 test_force_list();
1133 
1134 test_p_i_sim := procedure() {
1135   d := decision_algorithm.new();
1136   d.p := { ['A', 'P(A)'] };
1137   d.r := { ['P(A)', {0, 0.25, 0.75, 1}] };
1138   ds1 := { ['P(A)', 0.25] };
1139   ds2 := { ['P(A)', 0] };
1140   assert(abs(d.p_i_sim(ds1, ds2, 'P(A)') - 2 / 3) < 0.01,
1141          "Error: p_i_sim is incorrect.");
1142   print('[ok] p_i_sim');
1143 };
1144 test_p_i_sim();
1145 
1146 test_id_expr := procedure() {
1147   d := decision_algorithm.new();
1148   d.i := {'P(I)'};
1149   d.p := {['A', 'P(A)'], ['I', 'P(I)'] };
1150   assert(d.id_expr('A') == 'A' &&
1151          d.id_expr('I') == At_t('I', -1) &&
1152          d.id_expr(At_t('I', -1)) == At_t('I', -2),
1153          "Error: id_expr is incorrect.");
1154   print('[ok] id_expr');
1155 
1156   d.p := d.p + { [At_t('I', -1), 'P(I_-1)'],
1157                  [At_t('I', -2), 'P(I_-2)'] };
1158   assert(d.id_var('P(I)')    == 'P(I_-1)' &&
1159          d.id_var('P(I_-1)') == 'P(I_-2)',
1160          "Error: id_var is incorrect.");
1161   print('[ok] id_var');
1162 };
1163 test_id_expr();
1164 test_id_var := procedure() { test_id_expr(); };
1165 
1166 test_poss_i_events := procedure() {
1167   d := decision_algorithm.new();
1168   d.i := { 'P(I)', 'P(I2)' };
1169   d.r := { ['P(I)', {0,1}], ['P(I2)', {0,1}] };
1170   assert(d.poss_i_events() == {
1171            {["P(I)", 0]}, {["P(I)", 1]},
1172            {["P(I2)", 0]}, {["P(I2)", 1]},
1173            {["P(I)", 0], ["P(I2)", 0]},
1174            {["P(I)", 0], ["P(I2)", 1]},
1175            {["P(I)", 1], ["P(I2)", 0]},
1176            {["P(I)", 1], ["P(I2)", 1]}
1177          }, "Error: poss_i_events is incorrect.");
1178   print('[ok] poss_i_events');
1179 };
1180 test_poss_i_events();
1181 
1182 test_e_vars := procedure() {
1183   d := decision_algorithm.new();
1184   d.n := [{ ['e', {['A','E_n1(A)']}], ['m2', {'M2_n1'}], ['o', {'O_n1'}] },
1185           { ['e', {['B','E_n2(B)']}], ['m2', {'M2_n2'}], ['o', {'O_n2'}] } ];
1186   assert(d.e_vars() == { 'E_n1(A)', 'E_n2(B)' },
1187          "Error: e_vars  is incorrect.");
1188   print('[ok] e_vars ');
1189 };
1190 test_e_vars();
1191 
1192 test_poss_ue_events := procedure() {
1193   d := decision_algorithm.new();
1194   d.n := [{ ['e', {['A','E']}], ['m2', {'M2'}], ['o', {'O'}] }];
1195   d.u := { ['B', 'U'] };
1196   d.r := { ['E', {0,1}], ['U', {0,1}] };
1197   assert(d.poss_ue_events() == {
1198            {["E", 0]}, {["E", 0], ["U", 0]}, {["E", 0], ["U", 1]}, {["E", 1]},
1199            {["E", 1], ["U", 0]}, {["E", 1], ["U", 1]}, {["U", 0]}, {["U", 1]}
1200          }, "Error: poss_ue_events is incorrect.");
1201   print('[ok] poss_ue_events');
1202 };
1203 test_poss_ue_events();
1204 
1205 test_poss_p_i_events := procedure() {
1206   d := decision_algorithm.new();
1207   d.p := { ['A', 'P(A)'], ['B', 'P(B)'] };
1208   d.r := { ['P(A)', {0,1}], ['P(B)', {0,1}] };
1209   assert(d.poss_p_i_events() == { ["P(A)", 0], ["P(A)", 1], ["P(B)", 0], ["P(B)", 1] },
1210          "Error: poss_p_i_events is incorrect.");
1211   print('[ok] poss_p_i_events');
1212 };
1213 test_poss_p_i_events();
1214 
1215 test_poss_rx_states := procedure() {
1216   d := decision_algorithm.new();
1217   d.n := [{ ['o', {'O_n1'}] }, { ['o', {'O_n2'}] }];
1218   d.r := { ['O_n1', tf()], ['O_n2', tf()] };
1219   assert(d.poss_rx_states() == {
1220            {["O_n1", false], ["O_n2", false]},
1221            {["O_n1", false], ["O_n2",  true]},
1222            {["O_n1",  true], ["O_n2", false]},
1223            {["O_n1",  true], ["O_n2",  true]}
1224          }, "Error: poss_rx_states is incorrect.");
1225   print('[ok] poss_rx_states');
1226 };
1227 test_poss_rx_states();
1228 
1229 test_poss_rxs := procedure() {
1230   d := decision_algorithm.new();
1231   d.n := [{ ['o', {'O_n1_1', 'O_n1_2'}] }, { ['o', {'O_n2'}] } ];
1232   d.r := { [no_var, tf()] : no_var in {'O_n1_1', 'O_n1_2', 'O_n2'} };
1233   assert(d.poss_rxs() == {
1234            {["O_n1_1", false], ["O_n1_2", false]},
1235            {["O_n1_1", false], ["O_n1_2",  true]},
1236            {["O_n1_1",  true], ["O_n1_2", false]},
1237            {["O_n1_1",  true], ["O_n1_2",  true]},
1238            {["O_n2", false]},
1239            {["O_n2",  true]}
1240          }, "Error: poss_rxs is incorrect.");
1241   print('[ok] poss_rxs');
1242 };
1243 test_poss_rxs();
1244 
1245 /*
1246 test_n_disj := procedure() {
1247   d := decision_algorithm.new();
1248   b := causal_markov_model.new();
1249   b.u := { 'U1' };
1250   b.r := { ['U1', tf()] };
1251   be := { ['U1', true] };
1252   not_be := { ['U1', false] };
1253   w := b;
1254   d.cached_ref := { [be, { be, not_be }] };
1255   assert(d.n_disj(w, b, be) == 2,
1256          "Error: n_disj is incorrect.");
1257   print('[ok] n_disj');
1258 
1259   assert(d.n_ev_disj(w, b, not_be, be) == 1,
1260          "Error: n_ev_disj is incorrect.");
1261   print('[ok] n_ev_disj');
1262 
1263   assert(d.spec(w, b, not_be, be) == 1/2,
1264          "Error: spec is incorrect.");
1265   print('[ok] spec');
1266 };
1267 test_n_disj();
1268 test_n_ev_disj := procedure() { test_n_disj(); };
1269 test_spec := procedure() { test_spec(); };
1270 */
1271 
1272 test_ltr_up_to := procedure() {
1273   assert(da().ltr_up_to('A', 3) == { 'A_1', 'A_2', 'A_3' },
1274          "Error: ltr_up_to is incorrect.");
1275   print('[ok] ltr_up_to');
1276 };
1277 test_ltr_up_to();
1278 
1279 test_vars_up_to := procedure() {
1280   assert(da().vars_up_to({ 'A', 'B' }, 2) ==
1281            { 'A_1', 'A_2', 'B_1', 'B_2' },
1282          "Error: vars_up_to is incorrect.");
1283   print('[ok] vars_up_to');
1284 };
1285 test_vars_up_to();
1286 
1287 test_expr_to_ltr_var := procedure() {
1288   assert(da().expr_to_ltr_var(Implies('A', 'B'), 'U') == 'U(Implies(A, B))',
1289          "Error: expr_to_ltr_var is incorrect.");
1290   print('[ok] expr_to_ltr_var');
1291 };
1292 test_expr_to_ltr_var();
1293 
1294 test_poss_r_z  := procedure() {
1295   poss_z_vals := { { j : j in [1..x] } : x in [1..3] };
1296   assert(da().poss_r_z(poss_z_vals, da().ltr_up_to('A', 2)) == {
1297            {["A_1", {1}],       ["A_2", {1}]},
1298            {["A_1", {1}],       ["A_2", {1, 2}]},
1299            {["A_1", {1}],       ["A_2", {1, 2, 3}]},
1300            {["A_1", {1, 2}],    ["A_2", {1}]},
1301            {["A_1", {1, 2}],    ["A_2", {1, 2}]},
1302            {["A_1", {1, 2}],    ["A_2", {1, 2, 3}]},
1303            {["A_1", {1, 2, 3}], ["A_2", {1}]},
1304            {["A_1", {1, 2, 3}], ["A_2", {1, 2}]},
1305            {["A_1", {1, 2, 3}], ["A_2", {1, 2, 3}]}},
1306          "Error: poss_r_z is incorrect.");
1307   print('[ok] poss_r_z ');
1308 };
1309 test_poss_r_z();
1310 
1311 test_poss_r_ue := procedure() {
1312   u := { ['A', 'U(A)'], ['B', 'U(B)'] };
1313   poss_ue_vals := pow({ j : j in [1..2**2] }) - {{}};
1314   assert({ {["U(A)", {1, 2, 4}], ["U(B)", {1, 2, 4}]},
1315            {["U(A)", {2, 3, 4}], ["U(B)", {2, 3, 4}]}
1316          } < da().poss_r_ue(poss_ue_vals, u),
1317          "Error: poss_r_ue is incorrect.");
1318   print('[ok] poss_r_ue');
1319 };
1320 test_poss_r_ue();
  1 class causal_markov_model(n, u, v, r, f, a, cached_cm) {
  2   n := n; // integer number of time slices
  3   u := u; // Set of exogenous variables
  4           //   e.g. u := { 'O' };
  5   v := v; // Set of endogenous variables
  6           //   e.g. v := { 'F', 'L', 'ML' };
  7   r := r; // Function assigning a range of values for each variable
  8           //   e.g. r := { [x, { false, true }] : x in u + v };
  9   /* Higher order function from y in v to a higher order function from a
 10      function assigning values to exogenous and endogenous vars to a value 
 11      for y.
 12        eg f := { [y, f_y], [z, f_z], ... }
 13          f_y := { [f_y_1, true], [f_y_2, false], ... } 
 14           f_y_1 := { [y, false], [z, false], ... }       (row in truth table)
 15      A variable (technically, its past time slice) can also be its own parent
 16   */
 17   f := f;
 18 
 19   /* A list of actual states, one for each time slice. Each actual state is a 
 20      function assigning values to cm() vars of the given time.
 21        e.g. [ {['y_t1',true],...}, {['y_t2',false],...}, ... ] */
 22   a := a;
 23 
 24   uv := procedure() {
 25     return u + v;
 26   };
 27 
 28   u_t := procedure(t) {
 29     return { x + '_t' + t : x in u };
 30   };
 31   v_t := procedure(t) {
 32     return { x + '_t' + t : x in v };
 33   };
 34   uv_t := procedure(t) {
 35     return u_t(t) + v_t(t);
 36   };
 37 
 38   cm := procedure() {
 39     if (cached_cm == om) {
 40       cached_cm := compute_cm();
 41     }
 42     return cached_cm;
 43   };
 44 
 45     compute_cm := procedure() {
 46       return causal_model(cm_u(), cm_v(), cm_r(), cm_f(), om);
 47     };
 48 
 49     cm_u := procedure() {
 50       u_ts := { u_t(t) : t in [1..n] };
 51       if (u_ts == {}) {
 52         u_ts := {{}};
 53       }
 54       return +/ u_ts + v_t(1);
 55       // Includes the first v_t since there's no previous time parents
 56     };
 57 
 58     cm_v := procedure() {
 59       return +/ { v_t(t) : t in [2..n] };
 60       // Starts at 2 because the first, parentless v_t is counted in cm_u
 61     };
 62 
 63     cm_r := procedure() {
 64       return { [x, r[orig_var(x)]] : x in cm_u() + cm_v() };
 65     };
 66 
 67     cm_f := procedure() {
 68       return { [x, fix_time(x, f[orig_var(x)], var_time(x) - 1)] : x in cm_v() };
 69     };
 70 
 71       fix_time := procedure(x, f_x, t) {
 72         // Add time subscript to f_x
 73         f_prev_t := { [assign_time(g, t), f_x[g]] : g in domain(f_x) };
 74         // Extend from one timeslice to states over all vars but x
 75         ext_xs := (cm_u() + cm_v()) - uv_t(t) - {x}; // extended vars
 76         poss_ext_states := possible_states(ext_xs, cm_r()); // extended states
 77         // Output remains f_prev_t[prev_t_st] since it's a markov model
 78         return +/ { { [prev_t_st + ext_state, f_prev_t[prev_t_st] ]
 79                       : ext_state in poss_ext_states }
 80                     : prev_t_st in domain(f_prev_t)
 81                   };
 82       };
 83 
 84   poss_states := procedure() {
 85     return possible_states(u + v, r);
 86   };
 87 
 88   poss_events := procedure() {
 89     return +/ { 2 ** ps : ps in poss_states() };
 90   };
 91 
 92   actual_events := procedure() {
 93     return pow(actual_state()) - {{}};
 94   };
 95 
 96   actual_state := procedure() {
 97     return +/ { actual_state : actual_state in a };
 98   };
 99 
100   /* Set of all actual synchronic events, 
101      i.e. those taking place at a single time. */
102   actual_sync_events := procedure() {
103     return +/ { pow(actual_state) : actual_state in a } - {{}};
104   };
105 
106   /* Use the causal model on the first two times */
107   // g needs to be full state?
108   response :=  procedure(ys, g) {
109     g1 := assign_time(g, 1);
110     ys2 := { y + '_t' + 2 : y in ys };
111     return drop_time(cm().response(ys2, g1));
112   };
113 
114   /* Form the full state transition table and convert it to a string along with
115      the vars and ranges. */
116   t_to_str := procedure() {
117     tt := { [ps, response(v, ps)] : ps in poss_states() };
118     return str([tt, n, u, v, r]);
119   };
120 
121   static {
122     new := procedure() {
123       return causal_markov_model(0, {}, {}, {}, {}, [], om);
124     };
125   }
126 }
  1 print("\nTesting Causal Markov Model");
  2 print(  "===========================");
  3 
  4 cmm1 := cachedProcedure() {
  5   cmm := causal_markov_model.new();
  6   cmm.n := 3;
  7   cmm.u := { 'U' };
  8   cmm.v := { 'V' };
  9   cmm.r := { ['U', tf()], ['V', tf()] };
 10   cmm.f := { // eg v_t2 = U_t1 xor V_t1
 11     ['V', { [{ ['U', false], ['V', false] }, false],
 12             [{ ['U', false], ['V',  true] },  true],
 13             [{ ['U',  true], ['V', false] },  true],
 14             [{ ['U',  true], ['V',  true] }, false]
 15           }]
 16   };
 17   cmm.a := {
 18     {['U_t1',  true], ['V_t1', false]},
 19     {['U_t2',  true], ['V_t2',  true]},
 20     {['U_t3', false], ['V_t3',  true]}
 21   };
 22   return cmm;
 23 };
 24 
 25 test_cm_u := procedure() {
 26   assert(cmm1().cm_u() == { 'U_t1', 'U_t2', 'U_t3', 'V_t1' },
 27          "Error: cm_u is incorrect.");
 28   print('[ok] cm_u');
 29 };
 30 test_cm_u();
 31 
 32 test_cm_v := procedure() {
 33   assert(cmm1().cm_v() == { 'V_t2', 'V_t3' },
 34          "Error: cm_v is incorrect.");
 35   print('[ok] cm_v');
 36 };
 37 test_cm_v();
 38 
 39 test_cm_r := procedure() {
 40   assert(cmm1().cm_r() == { ['U_t1', tf()], ['U_t2', tf()], ['U_t3', tf()],
 41                             ['V_t1', tf()], ['V_t2', tf()], ['V_t3', tf()] },
 42          "Error: is incorrect.");
 43   print('[ok] cm_r');
 44 };
 45 test_cm_r();
 46 
 47 test_cm_f := procedure() {
 48   assert( { [ { ["U_t1", false], ["U_t2", false], ["U_t3", false],
 49                 ["V_t1", false],                  ["V_t3", false]  }, false],
 50             [ { ["U_t1", false], ["U_t2", false], ["U_t3", false],
 51                 ["V_t1", false],                  ["V_t3",  true]  }, false]
 52           } < cmm1().cm_f()['V_t2'],
 53           "Error: cm_f is incorrect.");
 54   print('[ok] cm_f');
 55 };
 56 test_cm_f();
 57 
 58 test_poss_states := procedure() {
 59   assert(cmm1().poss_states() == {
 60            { ["U", false], ["V", false] },
 61            { ["U", false], ["V",  true] },
 62            { ["U",  true], ["V", false] },
 63            { ["U",  true], ["V",  true] }
 64          },
 65          "Error: poss_states is incorrect.");
 66   print('[ok] poss_states');
 67 };
 68 test_poss_states();
 69 
 70 test_poss_events := procedure() {
 71   assert(cmm1().poss_events() == { {},
 72            {["U", false]}, {["U", true]}, {["V", false]}, {["V", true]},
 73            {["U", false], ["V", false]},
 74            {["U", false], ["V",  true]},
 75            {["U",  true], ["V", false]},
 76            {["U",  true], ["V",  true]}
 77          },
 78          "Error: poss_events is incorrect.");
 79   print('[ok] poss_events');
 80 };
 81 test_poss_events();
 82 
 83 test_actual_events := procedure() {
 84   assert(// test a random event
 85          { ["U_t2", true], ["U_t3", false], ["V_t1", false], ["V_t2", true]
 86          } in cmm1().actual_events(),
 87          "Error: actual_events is incorrect.");
 88   print('[ok] actual_events');
 89   print('[ok] actual_state');
 90 };
 91 test_actual_events();
 92 test_actual_state := procedure() { test_actual_state(); };
 93 
 94 test_actual_sync_events := procedure() {
 95   assert(// test some random events
 96          !({ ["U_t2", true], ["U_t3", false], ["V_t1", false], ["V_t2", true]
 97            } in cmm1().actual_sync_events()) &&
 98          { {['U_t2', true]}, {['U_t3',false], ['V_t3',true]}
 99            } < cmm1().actual_sync_events() &&
100          !({['U_t1', false]} in cmm1().actual_sync_events()),
101          "Error: actual_sync_events is incorrect.");
102   print('[ok] actual_sync_events');
103 };
104 test_actual_sync_events();
105 
106 test_response := procedure() {
107   assert(cmm1().response({'V'}, {['U', true], ['V', false]})['V'] == true &&
108          cmm1().response({'V'}, {['U', true], ['V',  true]})['V'] == false,
109          "Error: response is incorrect.");
110   print('[ok] response');
111 };
112 test_response();
113 
114 test_t_to_str := procedure() {
115   assert(eval(cmm1().t_to_str()) == [
116            // tt
117            { [{["U", false], ["V", false]}, {["V", false]}],
118              [{["U", false], ["V",  true]}, {["V",  true]}],
119              [{["U",  true], ["V", false]}, {["V",  true]}],
120              [{["U",  true], ["V",  true]}, {["V", false]}] },
121            3, // n
122            {"U"},
123            {"V"},
124            {["U", {false, true}], ["V", {false, true}]} // r
125          ],
126          "Error: t_to_str is incorrect.");
127   print('[ok] t_to_str');
128 };
129 test_t_to_str();
  1 /* Adapted from Judea Pearl's Causal Models
  2      e.g. Pearl, Judea. Causes and Explanations: A Structural-Model Approach
  3           http://ftp.cs.ucla.edu/pub/stat_ser/R266-part1.pdf
  4 */
  5 class causal_model(u, v, r, f, cached_poss_states) {
  6   u := u; // Set of exogenous variables
  7           //   e.g. u := { 'O' };
  8   v := v; // Set of endogenous variables
  9           //   e.g. v := { 'F', 'L', 'ML' };
 10   r := r; // Function assigning a range of values for each variable
 11           //   e.g. r := { [x, { false, true }] : x in u + v };
 12   /* Higher order function from y in v to a higher order function from a
 13      function assigning values to exogenous and other endogenous vars to a 
 14      value for y.
 15        eg f := { [y, f_y], [z, f_z], ... }
 16             f_y := { [f_y_1, true], [f_y_2, false], ... } // rows in truth table
 17               f_y_1 := { [x, false], [z, false], ... }
 18   */
 19   f := f;
 20 
 21   uv := procedure() {
 22     return u + v;
 23   };
 24 
 25   /* Takes an assignment of values to endogenous variables and returns a 
 26      causal submodel. See 7.12 Pearl's Causality book.  */
 27   effect_of_do := procedure(g) {
 28     h := { [y, f[y]] : y in v | !(y in domain(g)) } +
 29          { [x, i_proc_to_func(x, procedure(h) {
 30                                    return g[x];
 31                                  })] : x in domain(g) };
 32     return causal_model(u, v, r, h, om);
 33   };
 34 
 35   /* Takes a set of variables ys and an assignment of values g to endogenous 
 36      and exogenous variables and returns an assignment of values to ys */
 37   response := procedure(ys, g) {
 38     sm := effect_of_do(g);
 39     result := { [y, sm.solve_for(y, g)] : y in ys };
 40     if (dev) { print(result); }
 41     return result;
 42   };
 43 
 44   /* Takes a variable y and an assignment of values g to exogenous variables 
 45      and returns a value for y */
 46   solve_for := procedure(y, g) {
 47     g := { [x, g[x]] : x in domain(g) | x in u };
 48     if (y in u) {
 49       return g[u];
 50     } else if (#parents(y) == 0) {
 51       return last(first(f[y]));
 52     } else { // y in v
 53       // todo: sometimes a conflict? 
 54       h := g + { procedure() {
 55                    if (z in parents(y)) {
 56                      return [z, solve_for(z, g)];
 57                    } else { // doesn't matter, pick first
 58                      return [z, first(r[z])];
 59                    }
 60                  }() : z in (u + v) - domain(g) | z != y };
 61       return f[y][h];
 62     }
 63   };
 64 
 65   // Takes a variable and returns its parents
 66   parents := procedure(y) {
 67     return { z : z in u + v | z != y && is_nec_parent_of(z, y) };
 68   };
 69 
 70     /* Check whether there exists an assignment of values to the other 
 71        variables such that there are still different values of y depending 
 72        on values of z.  */
 73     is_nec_parent_of := procedure(z, y) {
 74       xs := {x : x in u + v | x != y && x != z}; // all other variables but z
 75       g_xs := possible_states(xs, r);
 76       return exists(g_x in g_xs | values_diff(g_x, y, z));
 77     };
 78 
 79     values_diff := procedure(g_x, y, z) {
 80       vals := { f[y][ g_x + {[z,r_i]} ] : r_i in r[z] };
 81       return #vals > 1;
 82     };
 83 
 84   /* Instance method shorthands for call to static version */
 85   i_proc_to_func := procedure(y, p) {
 86     return causal_model.proc_to_func(y, u, v, r, p);
 87   };
 88 
 89   poss_states := procedure() {
 90     if (cached_poss_states == om) {
 91       cached_poss_states := compute_poss_states();
 92     }
 93     return cached_poss_states;
 94   };
 95 
 96     compute_poss_states := procedure() {
 97       return possible_states(u + v, r);
 98     };
 99 
100   poss_events := cachedProcedure() {
101     return +/ { pow(ps) : ps in poss_states() };
102   };
103 
104   sorted_v := procedure() {
105     return sort_list( [y : y in v],
106                       procedure(y, z) { return depth(y) < depth(z); } );
107   };
108 
109     depth := procedure(y) {
110       if (y in u) {
111         return 0;
112       } else {
113         return 1 + max({ depth(z) : z in parents(y) });
114       }
115     };
116 
117   static {
118     /* Helper method that takes an endogenous variable y, & u, v, r and a 
119        procedure and returns a function. The procedure should take an 
120        assignment of values to exogenous variables and other endogenous 
121        variables and return y's value. */
122     proc_to_func := procedure(y, u, v, r, p) {
123       dom := possible_states({ x : x in u + v | x != y }, r);
124       return { [g, p(g)] : g in dom };
125     };
126 
127     new := procedure() {
128       return causal_model({}, {}, {}, {}, om);
129     };
130   }
131 }
 1 print("\nTesting Causal Model");
 2 print(  "====================");
 3 
 4 cm1 := procedure() {
 5   u := { 'O' };
 6   v := { 'F', 'L', 'ML' };
 7   r := { [x, { false, true }] : x in u + v };
 8   f_f := causal_model.proc_to_func('F', u, v, r, procedure(h) {
 9            if (h['L'] || h['ML']) {
10              return true;
11            } else {
12              return false;
13            }
14          });
15   f := { ['F', { [{ ['O', false], ['L', false], ['ML', false] }, false],
16                  [{ ['O', false], ['L', false], ['ML', true ] }, true ],
17                  [{ ['O', false], ['L', true ], ['ML', false] }, true ],
18                  [{ ['O', false], ['L', true ], ['ML', true ] }, true ],
19                  [{ ['O', true ], ['L', false], ['ML', false] }, false],
20                  [{ ['O', true ], ['L', false], ['ML', true ] }, true ],
21                  [{ ['O', true ], ['L', true ], ['ML', false] }, true ],
22                  [{ ['O', true ], ['L', true ], ['ML', true ] }, true ] }],
23          ['L',  causal_model.proc_to_func('L', u, v, r, procedure(g) {
24                   return false;
25                 })],
26          ['ML', causal_model.proc_to_func('ML', u, v, r, procedure(g) {
27                   return false;
28                 })] };
29   assert(f_f == f['F'], 'proc_to_func helper method is incorrect');
30   print('[ok] proc_to_func');
31 
32   cm := causal_model(u,v,r,f);
33   assert(cm.f['F'] == f_f, 'causal model class is incorrect');
34 
35   return cm;
36 };
37 test_proc_to_func := procedure() { cm1(); };
38 
39 test_effect_of_do := procedure() {
40   sm := cm1().effect_of_do({ ['F', true], ['L', true] });
41   assert((sm.f['F' ])[{ ['O', false], ['L', false], ['ML', false] }] <==> true,
42            'effect_of_do F is incorrect');
43   assert((sm.f['L' ])[{ ['O', false], ['F', false], ['ML', false] }] <==> true,
44            'effect_of_do L is incorrect');
45   assert((sm.f['ML' ])[{ ['O', false], ['F', false], ['L', false] }] <==> false,
46            'effect_of_do ML is incorrect');
47   print('[ok] effect_of_do');
48 };
49 
50 test_parents := procedure() {
51   assert(cm1().parents('F') == { 'L', 'ML' }, 'parents is incorrect');
52   print('[ok] parents');
53 };
54 
55 test_response := procedure() {
56   assert(cm1().response({'F'}, { ['L', false], ['O', true] })['F'] == false,
57          'response is incorrect');
58   assert(cm1().response({'F'}, { ['L', true], ['O', true] })['F'] == true,
59          'response is incorrect');
60   print('[ok] response');
61 };
  1 dev := false;
  2 
  3 tf := procedure() { return { true, false }; };
  4 
  5 // Non-empty powerset. Powerset without the empty set.
  6 ne_pow := procedure(s) {
  7   return pow(s) - {{}};
  8 };
  9 
 10 min := procedure(l) {
 11   mins := { e : e in l | forall(e2 in [x :  x in l | x != e ] | e2 >= e) };
 12   return first(mins);
 13 };
 14 
 15 max := procedure(l) {
 16   maxs := { e : e in l | forall(e2 in [x :  x in l | x != e ] | e >= e2) };
 17   return first(maxs);
 18 };
 19 
 20 avg := procedure(l) {
 21   return (+/ l)  / #l;
 22 };
 23 
 24 index_of := procedure(e, xs) {
 25   i := 1;
 26   while (xs[i] != e) {
 27     i := i + 1;
 28   }
 29   return i;
 30 };
 31 
 32 equiv_expr := procedure(e1, e2, dom) {
 33   if (dev) { print('equiv_expr dom ' + dom); }
 34   if (isList(e1) && !isList(e2) ||
 35       isList(e2) && !isList(e1)   ) {
 36     return false;
 37   } else if (isList(e1)) {
 38     return forall(j in [1..#e1] | equiv_expr(e1[j], e2[j], dom) );
 39   } else {
 40     // todo: handle quantifiers, []->?
 41     return expr_to_events(e1, dom) == expr_to_events(e2, dom);
 42   }
 43 };
 44 
 45 // Possible references for arity
 46 old_range_for_arity := cachedProcedure(e_arity, poss_events) {
 47   /* Singleton events, each being a set of just one event, a trivial 
 48      disjunction of truth-making events given our convention here of using 
 49      sets of events to express the disjunction of them. */
 50   sngl_events := { {event} : event in poss_events };
 51   if (e_arity == 0) {
 52     return sngl_events;
 53   /* Wrapping the sngl event below as a list/1-tuple might be a little 
 54      counterintuitive but it allows it to be treated consistently with
 55      higher arities */
 56   } else if (e_arity == 1) {
 57       return 2 ** { [se] : se in sngl_events };
 58   } else {
 59       return 2 ** iter_prod(sngl_events, e_arity);
 60   }
 61 };
 62 
 63 // Possible references for arity of expr
 64 range_for_arity := cachedProcedure(e_arity, poss_sts) {
 65   /* Propositions are sets of possible states. They can also be treated as 
 66      objects. You might think of them as the essential realization conditions
 67      for objects. */
 68   poss_objects := pow(poss_sts);
 69   // handling empty sets?
 70   if (e_arity == 0) {
 71     return poss_objects;
 72   /* Wrapping the sngl event below as a list/1-tuple might be a little 
 73      counterintuitive but it allows it to be treated consistently with
 74      higher arities */
 75   } else if (e_arity == 1) {
 76       return pow({ [x] : x in poss_objects - {{}} });
 77   } else {
 78       return pow(iter_prod(poss_objects - {{}}, e_arity));
 79       // return 2 ** iter_prod(sngl_events, e_arity);
 80   }
 81 };
 82 
 83 connectives := procedure() {
 84   return {'And','Not','Or','Implies','BoxArrow','At_t','Forall','Exists'};
 85 };
 86 
 87 /* Unused?
 88 fcts_n_strs := procedure(expr) {
 89   if (isTerm(expr)) {
 90     return ({ fct(expr) } - connectives()) + 
 91            (+/ { fcts_n_strs(expr) : arg_i in args(expr) });
 92   } else {
 93     return { expr };
 94   }
 95 };
 96 */
 97 
 98 maybe_not := procedure(var, val) {
 99   if (val == true) {
100     return var;
101   } else {
102     return Not(var);
103   }
104 };
105 
106 inv_maybe_not := procedure(expr) {
107   if (isTerm(expr) && fct(expr) == "Not") {
108     return { [expr, false] };
109   } else {
110     return { [expr,  true] };
111   }
112 };
113 
114 /* Takes an expr and returns the set of vars, atomic propositions and objects
115    that occur in it.  */
116 atomic_exprs := procedure(expr) {
117   if (isTerm(expr)) {
118     return +/ { atomic_exprs(arg_i) : arg_i in args(expr) };
119   } else {
120     return { expr };
121   }
122 };
123 
124 /* Takes an expr and returns the set of atomic propositions and objects
125    that occur in it (but not its variables).  */
126 atomic_terms := procedure(expr) {
127   return atomic_exprs(expr) - parse_vars(expr);
128 };
129 
130 /* Takes a function assigning true or false to vars and turns it into an 
131    expression that is a conjunction of either var or Not(var) depending 
132    on each variable's truth value.  */
133 func_to_expr := procedure(g) {
134   sorted_vs := sort([ v : v in domain(g) ]);
135   return iter_fct("And", [ maybe_not(v, g[v]) : v in sorted_vs ]);
136 };
137 
138 all_atomic_formulas := procedure(exprs) {
139   return +/ { atomic_formulas(expr) : expr in exprs };
140 };
141 
142 atomic_formulas := procedure(expr) {
143   return atomic_forms(expr) - parse_vars(expr);
144 };
145 
146 atomic_forms := procedure(expr) {
147   if (dev) { print('atomic forms ' + expr); }
148   if (!isTerm(expr) || !(fct(expr) in connectives())) {
149     return { expr };
150   } else {
151     return +/ { atomic_forms(arg_i) : arg_i in args(expr) };
152   }
153 };
154 
155 expr_to_events := procedure(expr, dom) {
156   state_space := exprs_to_state_space(dom);
157   if (dev) { print('st sp: ' + state_space); }
158   expand := procedure(events) {
159     if (dev) { print('expand ' + events); }
160     return { state : state in state_space
161              | exists(event in events | restricted_eq(state, event)) };
162   };
163   ete := procedure(expr) {
164     if (dev) { print('ete(' + expr + ')'); }
165     if (isString(expr)) {
166       return { state : state in state_space | state[expr] == true };
167     }
168     if (isTerm(expr)) {
169       if (#args(expr) >= 1) { a1 := args(expr)[1]; }
170       if (#args(expr) >= 2) { a2 := args(expr)[2]; }
171 
172       if (!(fct(expr) in connectives())) {  // applied predicate / relation
173         return { [expr, true] };
174       } else if (fct(expr) == "And") {
175         return expand(ete(a1)) * expand(ete(a2));
176         return expand(ete(a1)) * expand(ete(a2));
177         // return { g + h : [g,h] in ete(a1) >< ete(a2) | isMap(g + h) };
178       } else if (fct(expr) == "Or") {
179         return expand(ete(a1)) + expand(ete(a2));
180       // a -> b  ==  ~a v b
181       } else if (fct(expr) == "Implies") {
182         return expand(ete( Or(Not(a1), a2) ));
183       } else if (fct(expr) == "Exists") {
184         // Ex Qx <-> ~ Vx ~Qx
185         return expand(ete( Not( Forall(a1, Not(a2)) ) ));
186       } else if (fct(expr) == "Forall") {
187         terms := { var : var in dom | isString(var) };
188         return { state : state in state_space |
189           forall(term in terms | is_true_in(sub(a2, a1, term), state, dom)) };
190       } else if (fct(expr) == "Not") {
191         if (dev) { print('Not'); }
192         switch {
193           case isString(a1) :
194             if (dev) { print('Not(str)'); }
195             return { {[a1, false]} };
196           // Pair, i.e. { [var, val] }
197           case isSet(a1) && #a1 == 1 && isMap(a1) :
198             return { { [first(a1)[1], !first(a1)[2]] } };
199           case isSet(a1) && isMap(a1) :
200             if (dev) { print('~(a ^ b) = ~a v ~b'); }
201             // ~(a ^ b) = ~a v ~b
202             return +/ { expand(ete(Not({pair}))) : pair in a1 };
203           // Singleton map, eg { g } or { {[x,y],[w,z],...} }
204           case isSet(a1) && #a1 == 1 :
205             if (dev) { print('~ { g } = ~g'); }
206             // ~ { g } = ~g
207             return expand(ete(Not(first(a1))));
208           // Implicit disjunction of maps, eg { g, h, ... }
209           case isSet(a1) :
210             if (dev) { print('~(a v b) = ~a ^ ~b'); }
211             // ~(a v b) = ~a ^ ~b
212             head := first(a1);
213             return expand(ete(  And( Not(head), Not(a1 - {head}) )  ));
214           // Do we allow Not(a[]->b)? How?
215           // case isTerm(a1) && fct(a1) == 'BoxArrow' :
216           default : // e.g. functor term, i.e. predicate or other connective?
217             if (dev) { print('[~a] = [~[a]]'); }
218             // [~a] = [~[a]]
219             return expand(ete(Not( expand(ete(a1)))));
220         }
221       }
222     } else { // not term, atomic expr
223       return ({ [expr, true] });
224     }
225   }; // end ete
226   return expand(ete(expr));
227 };
228 
229 /* https://en.wikipedia.org/wiki/Fold_(higher-order_function)  */
230 foldr := procedure(proc, z, l) {
231   if (#l == 0) {
232     return z;
233   } else if (#l == 1) {
234     return proc(l[1], z);
235   } else {
236     return proc(l[1], foldr(proc, z, l[2..]));
237   }
238 };
239 
240 /* Takes a setlx functor and returns a procedure which takes two expressions
241    and if either is om, returns the other, or else constructs a term out of
242    the functor and two expressions.  */
243 fctize := procedure(fct) {
244   return procedure(e1, e2) {
245     if (e1 == om) { return e2; } else if (e2 == om) { return e1; } else {
246       return makeTerm(fct, [e1,e2]);
247     }
248   };
249 };
250 
251 /* Takes a setlx functor and list of expressions and constructs an expression
252    by repeatedly applying the functor.  */
253 iter_fct := procedure(fct, l) {
254   return foldr(fctize(fct), om, l);
255 };
256 
257 /* Takes a function f and returns its inverse. To avoid having to return a set 
258    or list, it just returns the first preimage that f maps to the img. */
259 inv := procedure(f) {
260   return { [img, ([ pre : pre in domain(f) | f[pre] == img])[1] ] : img in range(f) };
261 };
262 
263 /* Takes a setlx functor and returns its arity */
264 arity := procedure(fct) {
265   if (isTerm(fct)) {
266     return #args(fct);
267   } else {
268     return 0;
269   }
270 };
271 
272 /* Takes a set like s >< s >< s and returns tuples instead of nested pairs */
273 flatten := procedure(set) {
274   return { flatten_list(e) : e in set | isList(e)} +
275          { e : e in set | !isList(e) };
276 };
277 
278 flatten_list := procedure(list) {
279   if (isList(list)) {
280     return +/ [ flatten_list(l) : l in list ];
281   } else {
282     return [list];
283   }
284 };
285 
286 /* Takes a list of sets and returns the set of tuples taking a value from 
287    each set.  */
288 cart_prod := procedure(list) {
289   if (#list == 1) {
290     return first(list);
291   } else {
292     return flatten(first(list) >< cart_prod(list[2..]));
293   }
294 };
295 
296 /* Iterated Cartesian Product 
297    Yields the cartesian product of a set with itself n times.  */
298 iter_prod := procedure(set, n) {
299   if (n == 1) {
300     return set;
301   } else {
302     return flatten(iter_prod(set, n - 1) >< set);
303   }
304 };
305 
306 iter_mult := procedure(set) {
307   if (#set == 1) {
308     return first(set);
309   } else {
310     return first(set) * iter_mult(set - { first(set) });
311   }
312 };
313 
314 /* Takes a set of sets and returns a set of all possible choices taking an 
315    element from each set. */
316 choices := procedure(set) {
317   list := [ [elem : elem in subset] : subset in set ];
318   return set_choices(list);
319 };
320 
321 // Same as list_choices but returns a set rather than list of the choices.
322 set_choices := procedure(r) {
323   return { { o : o in p } : p in list_choices(r) };
324 };
325 
326 /* Takes a list of lists and returns a list of all possible choices taking
327    an element from each list  */
328 list_choices := procedure(r) {
329   if (#r == 1) {
330     return [[x] : x in r[1]];
331   } else {
332     return +/ [  [[x] + row : row in list_choices(r[2..])]  : x in r[1] ];
333   }
334 };
335 
336 /* Generate the set of all functions from dom to rng. 
337    Generate atomic maps. Collect the choices of those atomic maps. */
338 function_space := procedure(dom, rng) {
339   // Collect into a set the elements from a list of choices of atomic maps
340   // p for preimage and permutation? o is ordered pair?
341   return { { o : o in p } : p in list_choices([ [ [p, i] : i in rng ] : p in dom]) };
342 };
343 
344 /* Takes a set of variables s and a function r assigning a set which is the 
345    range of possible values to each variable and returns the set of possible 
346    functions assigning values to those variables.  */
347 possible_states := cachedProcedure(s, r) {
348   if (s == {} || r == {}) {
349     return {{}};
350   }
351   return choices({ {v} >< r[v] : v in s });
352 };
353 
354 /* Takes a function and returns all partial functions, including empty.  */
355 partial_functions := cachedProcedure(f) {
356   return 2 ** f;
357 };
358 
359 /* Takes two functions assigning values to vars and tests whether they are
360    compatible, i.e. do not assign different values to the same var.  */
361 compatible_values := procedure(f, g) {
362   h := f + g;
363   return forall(var in domain(f) + domain(g) | h[var] != om);
364 };
365 compat := procedure(f, g) { return compatible_values(f, g); }; // alias
366 
367 /* Takes partial state ps, which is a function assigning values to a set of
368    variables, and set of complete states cs, a set of functions assigning 
369    values to a larger set of variables and returns the subset of cs compatible 
370    with ps.  */
371 compatible_complete_states := cachedProcedure(ps, cs) {
372   return { c : c in cs
373            | forall(d in domain(ps)
374                     | c[d] == ps[d] )
375          };
376 };
377 compat_compl_sts := procedure(ps, cs) { // alias 
378   return compatible_complete_states(ps, cs);
379 };
380 
381 proc2func := procedure(p, dom) {
382   return { [d, p(d)] : d in dom };
383 };
384 
385 // Breaks apart a var label of form 'V_ti' into ['V', i]
386 parse_var := procedure(label) {
387   j := #label;
388   while (int(label[j]) != om && j != 1) {
389     j -= 1;
390   }
391   t := int(label[(j + 1)..]);
392   assert(t != om, "Error: Expected var label to have time subscript.");
393   var := label[1..(j - 2)];
394   return [var, t];
395 };
396 
397 orig_var := procedure(label) {
398   return parse_var(label)[1];
399 };
400 
401 var_time := procedure(label) {
402   return parse_var(label)[2];
403 };
404 
405 /* Takes an assignment of values to variables and adds a time subscript to the 
406    variables' labels. */
407 assign_time := procedure(g, t) {
408   return { [ y+'_t'+t, g[y] ] : y in domain(g) };
409 };
410 
411 /* Takes an assignment of values to variables and removes time subscripts 
412    from the variables' labels. */
413 drop_time := procedure(g) {
414   return { [ orig_var(y), g[y] ] : y in domain(g) };
415 };
416 
417 /* From Hilbert System, Wikipedia 
418      https://en.wikipedia.org/w/index.php?title=Hilbert_system&oldid=873786806
419 */
420 is_logical_truth := procedure(expr) {
421   vars := parse_vars(expr);
422 
423   /* Axiom schemes
424      Redundant: 
425        P1: phi -> phi
426        as_p1 := procedure(phi) { return Implies(phi, phi); };  
427       
428      P2: phi -> (psi -> phi)  */
429   as_p2 := procedure(phi, psi) { return Implies(phi, Implies(psi, phi)); };
430   // P3: (phi -> (psi -> xi)) -> ((phi -> psi) -> (phi -> xi))
431   as_p3 := procedure(phi, psi, xi) {
432     return Implies(
433              Implies(phi, Implies(psi, xi)),
434              Implies(Implies(phi, psi), Implies(phi, xi))
435            );
436   };
437   // P4: (~phi -> ~psi) -> (psi -> phi)
438   as_p4 := procedure(phi, psi) {
439     return Implies(Implies(Not(phi), Not(psi)), Implies(psi, phi));
440   };
441   // Q5: Vx(phi) -> phi[x := t] where t may be substituted for x in phi
442   as_q5 := procedure(var, phi) {
443     return { Implies( Forall(var, phi), sub(phi, var, term))
444              : term in atomic_terms(expr) };
445   };
446   // Q6: Vx(phi -> psi) -> (Vx(phi) -> Vx(psi))
447   as_q6 := procedure(var, phi, psi) {
448     return Implies( Forall(var, Implies(phi, psi)),
449                     Implies(Forall(var, phi), Forall(var, psi)) );
450   };
451   // Q7: phi -> Vx phi where x not in phi
452   as_q7 := procedure(phi, var) {
453     if (var in atomic_exprs(phi)) {
454       return {};
455     } else {
456       return { Implies(phi, Forall(var, phi)) };
457     }
458   };
459   /* Conservative Extensions
460      Conjunction:
461        Introduction: phi -> (psi -> (phi ^ psi))  */
462   conj_intr := procedure(phi, psi) {
463     return Implies(phi, Implies(psi, And(phi, psi)));
464   };
465   //   Elimination: (phi ^ psi) -> phi   and   (phi ^ psi) -> psi
466   conj_elim := procedure(phi, psi) {
467     return { Implies(And(phi, psi), phi),
468              Implies(And(phi, psi), psi) };
469   };
470   /* Disjunction:
471        Introduction: phi -> (phi v psi)   and   psi -> (phi v psi)  */
472   disj_intr := procedure(phi, psi) {
473     return { Implies(phi, Or(phi, psi)),
474              Implies(psi, Or(phi, psi)) };
475   };
476   //   Elimination: (phi -> xi) -> ((psi -> xi) -> ((phi v psi) -> xi))
477   disj_elim := procedure(phi, psi, xi) {
478     return Implies( Implies(phi, xi),
479                     Implies( Implies(psi, xi),
480                              Implies( Or(phi, psi), xi)));
481   };
482   /* Existential Quantification 
483      Introduction: Vx(phi -> Ey(phi[x := y]))  */
484   ex_intr := procedure(phi, x_var, y_var) {
485     return Forall(x_var, Implies(phi, Exists(y_var, sub(phi, x_var, y_var))));
486   };
487   // Elimination: Vx(phi -> psi) -> (Ex(phi) -> psi) where !(x in psi)
488   ex_elim := procedure(phi, psi, var) {
489     if (var in atomic_exprs(psi)) {
490       return {};
491     } else {
492       return { Implies( Forall(var, Implies(phi, psi)),
493                         Implies( Exists(var, phi), psi)) };
494     }
495   };
496   // Modus Ponens Rule of Inference
497   // phi, phi -> psi |- psi
498   modus_ponens := procedure(phi, psi) {
499     if (isTerm(psi) && fct(psi) == 'Implies' && args(psi)[1] == phi) {
500       return { args(psi)[2] };
501     } else {
502       return {};
503     }
504   };
505 
506   derive_more := procedure(exprs) {
507     return
508       (+/ { as_q7(phi, var) : [phi, var] in exprs >< vars }) +
509       (+/ { { as_p2(phi, psi), as_p4(phi, psi), conj_intr(phi, psi) } +
510               { as_q6(var, phi, psi) : var in vars } +
511               conj_elim(phi, psi) + disj_intr(phi, psi) + modus_ponens(phi, psi)
512             : [phi, psi] in exprs >< exprs }) +
513       { as_p3(phi, psi, xi) : [phi, psi, xi] in iter_prod(exprs, 3) } +
514       (+/ { as_q5(var, phi, vars) : [var, phi] in vars >< exprs }) +
515       (+/ { ex_elim(phi, psi, var)
516             : [phi, psi, var] in cart_prod([exprs, exprs, vars]) } ) +
517       { ex_intr(phi, x_var, y_var)
518             : [phi, x_var, y_var] in cart_prod([exprs,vars,vars]) };
519   };
520   derive := procedure(exprs, j) {
521     if (j == 0) {
522       return exprs;
523     } else {
524       return derive_more(derive(exprs, j - 1));
525     }
526   };
527   return expr in derive(subformulas(expr), 10 ** 6);
528 };
529 
530 parse_vars := procedure(expr) {
531   if (isTerm(expr)) {
532     parse_more := +/ { parse_vars(arg) : arg in args(expr) };
533     if (fct(expr) in {'Forall', 'Exists'}) {
534       return { args(expr)[1] } + parse_more;
535     } else {
536       return parse_more;
537     }
538   } else {
539     return {};
540   }
541 };
542 
543 /* Substitute any occurrences of var in expr with sub_var instead. */
544 sub := procedure(expr, var, sub_var) {
545   if (isTerm(expr)) {
546     return makeTerm(fct(expr), [sub(arg, var, sub_var) : arg in args(expr)]);
547   } else {
548     if (expr == var) {
549       return sub_var;
550     } else {
551       return expr;
552     }
553   }
554 };
555 
556 subformulas := procedure(expr) {
557   if (isTerm(expr) && fct(expr) in connectives()) {
558     return { expr } + (+/ { subformulas(arg) : arg in args(expr) })
559                     - parse_vars(expr);
560   } else {
561     return { expr };
562   }
563 };
564 
565 /* Hausdorff Distance
566    "It is the greatest of all the distances from a point in one set to the 
567     closest point in the other set."   
568       -- https://en.wikipedia.org/wiki/Hausdorff_distance 
569 */
570 hausdorff_dist := procedure(s1, s2, dist_f) {
571   dist_to_closest := procedure(e1) {
572     return min([ dist_f(e1, e2) : e2 in s2 ]);
573   };
574   return max([ dist_to_closest(e1) : e1 in s1 ]);
575 };
576 
577 /* Kendalls Tau Distance
578    Takes two strict ordinal utility functions with the same domain and counts 
579    the number of discordant pairs. Returns a normalized distance from 0 - 1.
580    Each u1, u2 is a set of setlx functors Succ(state1, state2). (Succ stands 
581    for succeeds.) States are logical expressions, i.e. setlx functors of 
582    combinators on string base propositions. 
583 
584    "Kendall tau distance is also called bubble-sort distance since it is 
585     equivalent to the number of swaps that the bubble sort algorithm would 
586     take to place one list in the same order as the other list."
587       -- https://en.wikipedia.org/wiki/Kendall_tau_distance  
588 */
589 kendalls_tau_dist := procedure(u1, u2) {
590   pairs := { [x,y] : [x,y] in states_u(u1) >< states_u(u1) | x != y };
591   // Discordant (unordered) pairs
592   disc_un_pairs := #{ [x,y] : [x,y] in pairs
593                       | (Succ(x,y) in u1 && Succ(y,x) in u2) ||
594                         (Succ(x,y) in u2 && Succ(y,x) in u1)
595                     } / 2;
596   n := #states_u(u1);
597   return disc_un_pairs / (n * (n - 1) / 2);
598 };
599 
600 /* Takes an ordinal utility function expressed as a set of SuccEq setlx
601    functors comparing events or states and returns the set of atomic variables
602    those events or states range over.  */
603 atom_u := procedure(u) {
604   states := +/ { { args(pref)[1], args(pref)[2] } : pref in u };
605   return +/ { domain(state) : state in states };
606 };
607 
608 /* Takes an ordinal utility function and returns the states that are compared
609    in them.  */
610 states_u := procedure(u) {
611   return +/ { { args(pref)[1], args(pref)[2] } : pref in u };
612 };
613 
614 /* Cardinal Utility Function to Ordinal Utility Function */
615 card_u_to_ord_u := procedure(card_u, state_space) {
616   compare_sts := procedure(st1, st2) {
617     if (card_u[st1] >= card_u[st2]) {
618       return SuccEq(st1, st2); // Succeeds or equals
619     } else {
620       return SuccEq(st2, st1);
621     }
622   };
623   return { compare_sts(st1, st2) : [st1,st2] in state_space >< state_space };
624 };
625 
626 /* Additive Utility Function to Cardinal Utility Function: 
627    Takes a utility function expressed as { [expr, val], ... } where the 
628    utility of a state is the sum of all its true exprs' val.
629    Returns a cardinal utility function assigning values to states rather than
630    expressions.
631 */
632 add_u_to_card_u := procedure(add_u) {
633   atoms := +/ { atomic_formulas(expr) : expr in domain(add_u) };
634   state_space := exprs_to_state_space(domain(add_u));
635   util := procedure(state) {
636     sum := +/ { val : [expr, val] in add_u | is_true_in(expr, state, atoms) };
637     if (sum == om) { sum := 0; }
638     return sum;
639   };
640   return { [state, util(state)] : state in state_space };
641 };
642 
643 /* Expressions to State Space:
644    Return the possible ways of assigning true or false to each expr in exprs. */
645 exprs_to_state_space := procedure(exprs) {
646   if (dev) { print('exprs_to_state_space ' + exprs); }
647   atoms := +/ { atomic_formulas(expr) : expr in exprs };
648   return choices({ {[atom,true], [atom,false]} : atom in atoms });
649 };
650 
651 /* Expression is True in State:
652    Takes an expression and a state, a function assigning true/false to the 
653    atomic terms in expr, and returns true or false for whether the expr
654    is true in the state.
655 */
656 is_true_in := procedure(expr, state, dom) {
657   // return state in expr_to_events(expr, dom); // assumes event is state
658   return exists(event in expr_to_events(expr, dom)
659                 | restricted_eq(state, event) );
660 };
661 
662 /* Ordinality Utility Function Distance:
663    Takes two weak ordinal utility functions expressed as a set of SuccEq setlx 
664    functors and returns the distance between them.
665 
666    See hausdorff_dist and kendalls_tau_dist.
667    See also: https://en.wikipedia.org/wiki/Kemeny-Young_method 
668      Interestly, since we read the utility functions off the brains rather
669      than elicit them through votes, we don't need to be concerned with 
670      strategic manipulation through insincere votes.
671 */
672 ord_u_dist := procedure(u1, u2) {
673   atom_us := atom_u(u1) + atom_u(u2);
674   state_space := exprs_to_state_space(atom_us);
675   sharpenings := procedure(u) {
676     poss_total_orders := permutations([s : s in state_space]);
677     list_to_ord := procedure(l) {
678       if (l == []) {
679         return {};
680       } else {
681         return { Succ(l[1], s) : s in l[2..] } + list_to_ord(l[2..]);
682       }
683     };
684     ord_us := { list_to_ord(pto) : pto in poss_total_orders };
685     ext_u := atom_us - atom_u(u); // extending u vars
686     if (ext_u == {}) {
687       // return { u };
688       exts := {{}};
689     } else {
690       exts := choices({ {[atom,true], [atom,false]} : atom in ext_u });
691     }
692     return {
693       ord_u : ord_u in ord_us
694       | forall(pref in u // u: a1 >= a2 
695           | Succ(args(pref)[1], args(pref)[2]) in ord_u || // ord_u: a1 > a2
696             SuccEq(args(pref)[2], args(pref)[1]) in u || // u: a2 >= a1, ie ==
697             forall(ext in exts
698               | Succ(args(pref)[1] + ext, args(pref)[2] + ext) in ord_u
699                   // a1 + z > a2 + z
700             ) // end forall ext
701         ) // end forall pref
702     }; // end return
703   }; // end sharpenings  
704   return hausdorff_dist( sharpenings(u1),
705                          sharpenings(u2),
706                          kendalls_tau_dist );
707 };
708 
709 /* p 56 of setlx tutorial */
710 sort_list := procedure(l, cmp) {
711   if (#l < 2) { return l; }
712   m := #l \ 2;
713   [l1, l2] := [l[.. m], l[m+1 ..]];
714   [s1, s2] := [sort_list(l1, cmp), sort_list(l2, cmp)];
715   return merge(s1, s2, cmp);
716 };
717 
718 sort_set := procedure(s, cmp) {
719   return sort_list([x : x in s], cmp);
720 };
721 
722 merge := procedure(l1, l2, cmp) {
723   match ([l1, l2]) {
724     case [[], _] : return l2;
725     case [_, []] : return l1;
726     case [[x1|r1], [x2|r2]] :
727     if (cmp(x1, x2)) {
728       return [x1 | merge(r1, l2, cmp)];
729     } else {
730       return [x2 | merge(l1, r2, cmp)];
731     }
732   }
733 };
734 
735 setlx_chars := procedure() {
736   str := " !\"#%&'()*+,-./0123456789:;<=>?" +
737          "ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]_" +
738          "abcdefghijklmnopqrstuvwxyz{|}~";
739   return { str[j] : j in [1..#str] };
740 };
741 
742 all_strs_lte := procedure(s) {
743   poss_strs := {''};
744   length := 0;
745   while (length < #s) {
746     length := length + 1;
747     poss_strs := poss_strs + { poss_str + char : [poss_str, char] in
748                                                  poss_strs >< setlx_chars() };
749   }
750   return poss_strs;
751 };
752 
753 restrict_domain := procedure(event, vars) {
754   return { [var, val] : [var, val] in event | var in vars };
755 };
756 
757 restricted_eq := procedure(state, event) {
758   return event == restrict_domain(state, domain(event));
759 };
760 
761 /* Kolmogorov Complexity of a string
762    https://en.wikipedia.org/wiki/Kolmogorov_complexity  
763 
764    k is an incomputable function but we include a naive attempt to define it
765    anyway. Although it may seem possible this way, it will not work because
766    eval(poss_str) may not terminate. In practice, some finite approximation
767    could be used. For instance, if setlx allowed us to catch a timeout error
768    as many other languages do, it would be particularly easy by just trying 
769    each program string up to a fixed length of time.
770 */
771 k := procedure(s) {
772   poss_strs := {''};
773   success := false;
774   length := 0;
775   while (!success) {
776     length := length + 1;
777     poss_strs := { poss_str + char
778                    : [poss_str, char] in poss_strs >< setlx_chars() };
779     if ({ eval(poss_str) == s : poss_str in poss_strs } ||
780         length == #s + 2 /* '"'+s+'"' */                  )  {
781       success := true;
782     }
783   }
784   return length;
785 };
786 
787 /* Conditional Kolmogorov Complexity  K(s|i)  
788    The length of the shortest program which outputs a string s when given the 
789    string i as input. Again, this is incomputable but can be finitely 
790    approximated.
791 */
792 k_given := procedure(s, i) {
793   eval_str := procedure(str) {
794     try {
795       return eval(str + '("' + i + '");');
796     } catch(e) {
797       return '';
798     }
799   };
800   strs := { str : str in all_str_lte(#s + 2) | eval_str(str) == s };
801   strs := sort(strs, procedure(s1, s2) { return #s1 < #s2; });
802   return #strs[1];
803 };
  1 print('Testing Lib');
  2 print('===========');
  3 
  4 test_ne_pow := procedure() {
  5   assert(ne_pow({1,2}) == {{1}, {2}, {1,2}},
  6          "Error: ne_pow is incorrect.");
  7   print('[ok] ne_pow');
  8 };
  9 test_ne_pow();
 10 
 11 test_min := procedure() {
 12   assert(min([3,4,5]) == 3, "Error: min is incorrect.");
 13   print('[ok] min');
 14 };
 15 test_min();
 16 
 17 test_max := procedure() {
 18   assert(max([3,4,5]) == 5, "Error: max is incorrect.");
 19   print('[ok] max');
 20 };
 21 test_max();
 22 
 23 test_avg := procedure() {
 24   assert(avg([3,4,5]) == 4, "Error: avg is incorrect.");
 25   print('[ok] avg');
 26 };
 27 test_avg();
 28 
 29 test_index_of := procedure() {
 30   assert(index_of(5, [3,4,5]) == 3, "Error: index_of is incorrect.");
 31   print('[ok] index_of');
 32 };
 33 test_index_of();
 34 
 35 test_equiv_expr := procedure() {
 36   e1 := And('A', 'B');
 37   e2 := And('B', 'A');
 38   assert(equiv_expr(e1, e2, { 'A', 'B' }), "Error: equiv_expr is incorrect.");
 39   print('[ok] equiv_expr: A ^ B == B ^ A');
 40 
 41   assert(equiv_expr(Exists('x', Odd('x')),
 42                     Not(Forall('y', Not(Odd('y')))),
 43                     { Odd('a'), Odd('b'), Odd('c') }),
 44          "Error: equiv_expr is incorrect.");
 45   print('[ok] equiv_expr: Ex Odd(x) == ~( Vy ~Odd(y) )');
 46   // print('[ok] expr_to_events');
 47 };
 48 test_equiv_expr();
 49 
 50 test_maybe_not := procedure() {
 51   assert(maybe_not('A', true) == 'A',
 52          "Error: maybe_not true is incorrect.");
 53   print('[ok] maybe_not true');
 54   assert(maybe_not('A', false) == Not('A'),
 55          "Error: maybe_not false is incorrect.");
 56   print('[ok] maybe_not false');
 57 };
 58 test_maybe_not();
 59 
 60 test_atomic_exprs := procedure() {
 61   assert(atomic_exprs(Exists('x', And('B', Odd('a')))) == { 'x', 'B', 'a' },
 62          "Error: atomic_exprs is incorrect.");
 63   print('[ok] atomic_exprs');
 64 };
 65 test_atomic_exprs();
 66 
 67 test_atomic_terms := procedure() {
 68   assert(atomic_terms(Exists('x', And('B', Odd('a')))) == { 'B', 'a' },
 69          "Error: atomic_terms is incorrect.");
 70   print('[ok] atomic_terms');
 71 };
 72 test_atomic_terms();
 73 
 74 test_func_to_expr := procedure() {
 75   assert(func_to_expr({ ['A', true], ['B', false], ['C', false] }) ==
 76          And("A", And(Not("B"), Not("C"))),
 77          "Error: func_to_expr is incorrect.");
 78   print('[ok] func_to_expr');
 79   print('[ok] iter_fct');
 80   print('[ok] foldr');
 81   print('[ok] fctize');
 82 };
 83 test_func_to_expr();
 84 
 85 test_atomic_formulas := procedure() {
 86   assert(atomic_formulas(Forall('x', Implies(Odd('a'),'B') )) ==
 87          { Odd('a'), 'B' },
 88          "Error: atomic_formulas is incorrect.");
 89   print('[ok] atomic_formulas');
 90   print('[ok] atomic_forms');
 91 };
 92 test_atomic_formulas();
 93 test_atomic_forms := procedure() { test_atomic_formulas(); };
 94 
 95 test_expr_to_events := procedure() {
 96   assert(expr_to_events(Not('A'), {'A'}) == { { ['A', false] } },
 97          "Error: expr_to_events ~A is incorrect.");
 98   assert(expr_to_events('A', {'A', 'B' }) ==
 99          { { ['A', true], ['B', false] },
100            { ['A', true], ['B', true ] } },
101          "Error: expr_to_events A is incorrect.");
102   assert(expr_to_events(Implies('A', 'B'), {'A', 'B' }) ==
103          { { ['A', true ], ['B', true ] },
104            { ['A', false], ['B', true ] },
105            { ['A', false], ['B', false] } },
106          "Error: expr_to_events A->B is incorrect.");
107   dom := { 'A', 'B', 'C' };
108   assert(expr_to_events(And('A', 'B'), dom) ==
109          { { ['A', true], ['B', true], ['C', true ] },
110            { ['A', true], ['B', true], ['C', false] }  },
111          "Error: expr_to_events A^B is incorrect.");
112   print('[ok] expr_to_events');
113   print('[  ] more expr_to_events');
114 };
115 test_expr_to_events();
116 
117 test_inv := procedure() {
118   assert(inv({ ['a', 1], ['b', 2] }) == { [1, 'a'], [2, 'b'] },
119          "Error: inv is incorrect.");
120   print('[ok] inv');
121 };
122 test_inv();
123 
124 test_arity := procedure() {
125   assert(arity('A') == 0 && arity(B()) == 0 && arity(Odd('x')) == 1 &&
126          arity(GreaterThan('x', 'y')) == 2 &&
127          arity(Relation('x', 'y', 'z')) == 3,
128          "Error: arity is incorrect.");
129   print('[ok] arity');
130 };
131 test_arity();
132 
133 test_flatten := procedure() {
134   assert([1,'a','C'] in flatten({1, 2, 3} >< {'a','b','c'} >< {'A', 'B', 'C'}),
135          "Error: flatten is incorrect.");
136   print('[ok] flatten');
137   print('[ok] flatten_list');
138 };
139 test_flatten();
140 
141 test_cart_prod := procedure() {
142   assert([2,'a','B'] in cart_prod([{1, 2, 3}, {'a','b','c'}, {'A', 'B', 'C'}]),
143          "Error: cart_prod is incorrect.");
144   print('[ok] cart_prod');
145 };
146 test_cart_prod();
147 
148 test_iter_prod := procedure() {
149   assert([3,2,1] in iter_prod({1,2,3}, 3),
150          "Error: is incorrect.");
151   print('[ok] iter_prod ');
152 };
153 test_iter_prod();
154 
155 test_iter_mult := procedure() {
156   assert(iter_mult({ 3, 5, 7 }) == 105,
157          "Error: iter_mult is incorrect.");
158   print('[ok] iter_mult');
159 };
160 test_iter_mult();
161 
162 test_choices := procedure() {
163   assert({3,'b','A'} in choices({ {1,2,3}, {'a','b','c'}, {'A','B','C'} }),
164          "Error: is incorrect.");
165   print('[ok] choices ');
166   print('[ok] set_choices ');
167   print('[ok] list_choices ');
168 };
169 test_choices();
170 
171 test_function_space := procedure() {
172   dom := { 'a', 'b' };
173   rng := { 1, 2 };
174   assert(function_space(dom, rng) == { {["a", 1], ["b", 1]},
175                                        {["a", 1], ["b", 2]},
176                                        {["a", 2], ["b", 1]},
177                                        {["a", 2], ["b", 2]}
178                                      },
179          'Error: function_space is incorrect');
180   print('[ok] function_space');
181 };
182 test_function_space();
183 
184 test_possible_states := procedure() {
185   s := { 'a', 'b' };
186   r := { ['a', {0,1}], ['b', {3,4}] };
187   assert(possible_states(s,r) == { {["a", 0], ["b", 3]},
188                                    {["a", 0], ["b", 4]},
189                                    {["a", 1], ["b", 3]},
190                                    {["a", 1], ["b", 4]}
191                                  },
192          'Error: possible_states is incorrect');
193   print('[ok] possible_states');
194 };
195 test_possible_states();
196 
197 test_partial_functions := procedure() {
198   f := { [1,2], [2,4] };
199   assert(partial_functions(f) == {{}, {[1, 2]}, {[1, 2], [2, 4]}, {[2, 4]}},
200          'Error: partial_functions is incorrect');
201 };
202 test_partial_functions();
203 
204 test_compatible_values := procedure() {
205   assert(compatible_values(
206            { ['A',true], ['B',false] },
207            { ['C',false], ['B',false] }
208          ) == true,
209          "Error: compatible_values is incorrect.");
210   assert(compatible_values(
211            { ['A',true], ['B',false] },
212            { ['C',false], ['B',true] }
213          ) == false,
214          "Error: compatible_values is incorrect.");
215   print('[ok] compatible_values');
216 };
217 test_compatible_values();
218 
219 test_compatible_complete_states := procedure() {
220   ps := { ['X', true] };
221   cs := possible_states( { 'X', 'Y' },
222                          { ['X', {true, false}],
223                            ['Y', {true, false}]  });
224   assert(compatible_complete_states(ps, cs) == {
225          {["X", true], ["Y", false]},
226          {["X", true], ["Y", true]}
227        }, 'Error: compatible_complete_states is incorrect');
228   print('[ok] compatible_complete_states');
229 };
230 test_compatible_complete_states();
231 
232 test_proc2func := procedure() {
233   assert(proc2func(procedure(x) { return 2 * x; }, { 1,2,3 }) ==
234          { [1,2], [2,4], [3,6] },
235          "Error: is incorrect.");
236   print('[ok] proc2func');
237 };
238 test_proc2func();
239 
240 test_parse_var := procedure() {
241   assert(parse_var('V_t4') == ['V', 4],
242          "Error: parse_var is incorrect.");
243   print('[ok] parse_var');
244 };
245 test_parse_var();
246 
247 test_orig_var := procedure() {
248   assert(orig_var('V_t3') == 'V',
249          "Error: orig_var is incorrect.");
250   print('[ok] orig_var ');
251 };
252 test_orig_var();
253 
254 test_var_time := procedure() {
255   assert(var_time('V_t2') == 2,
256          "Error: is incorrect.");
257   print('[ok] var_time');
258 };
259 test_var_time();
260 
261 test_assign_time := procedure() {
262   assert(assign_time({ ['A',    true], ['B',     false] }, 4) ==
263                      { ['A_t4', true], ['B_t4', false] },
264          "Error: assign_time is incorrect.");
265   print('[ok] assign_time');
266 };
267 test_assign_time();
268 
269 test_drop_time := procedure() {
270   assert(drop_time({ ['A_t3', true], ['B_t3', false] }) ==
271                    { ['A',     true], ['B',    false] },
272          "Error: drop_time is incorrect.");
273   print('[ok] drop_time');
274 };
275 test_drop_time();
276 
277 /*
278 test_is_logical_truth := procedure() {
279   assert(true,
280          "Error: is_logical_truth is incorrect.");
281   print('[  ] is_logical_truth');
282 };
283 test_is_logical_truth();
284 */
285 
286 test_parse_vars := procedure() {
287   assert(parse_vars(Forall('x', Exists('y', And(Odd('b'), LessThan('y','a')))))
288          == { 'x', 'y' },
289          "Error: is incorrect.");
290   print('[ok] parse_vars');
291 };
292 test_parse_vars();
293 
294 test_sub := procedure() {
295   assert(sub(Forall('x', Exists('y', And(Odd('b'), LessThan('y','x')))),'x','z')
296          ==  Forall('z', Exists('y', And(Odd('b'), LessThan('y','z')))),
297          "Error: sub is incorrect.");
298   print('[ok] sub');
299 };
300 test_sub();
301 
302 test_subformulas := procedure() {
303   assert(subformulas(Forall('x', Exists('y', And(Odd('b'), LessThan('y','x')))))
304          == { Odd('b'), LessThan('y','x'), And(Odd('b'), LessThan('y','x')),
305               Exists('y', And(Odd('b'), LessThan('y','x'))),
306               Forall('x', Exists('y', And(Odd('b'), LessThan('y','x'))))
307             },
308          "Error: subformulas is incorrect.");
309   print('[ok] subformulas');
310 };
311 test_subformulas();
312 
313 test_hausdorff_dist := procedure() {
314   assert(hausdorff_dist({1, 8, 17}, {4, 7, 15},
315                         procedure(x, y) { return (x - y) ** 2; }) == 9,
316          "Error: hausdorff_dist is incorrect.");
317   print('[ok] hausdorff_dist');
318 };
319 test_hausdorff_dist();
320 
321 test_kendalls_tau_dist := procedure() {
322   // In actual usage, instead of letters, states/events would be compared.
323   // A > B > C > D
324   u1 := { Succ('A', 'B'), Succ('A', 'C'), Succ('A', 'D'),
325           Succ('B', 'C'), Succ('B', 'D'),
326           Succ('C', 'D') };
327   // C > A > D > B
328   // Bubble sort: abcd => aCBd => CAbd => caDB 
329   u2 := { Succ('C', 'A'), Succ('C', 'D'), Succ('C', 'B'),
330           Succ('A', 'D'), Succ('A', 'B'),
331           Succ('D', 'B') };
332   assert(kendalls_tau_dist(u1, u2) == 1/2, // = 3 / (4 * (4-1) / 2)
333          "Error: is incorrect.");
334   print('[ok] kendalls_tau_dist');
335   print('[ok] states_u');
336 };
337 test_kendalls_tau_dist();
338 
339 test_atom_u := procedure() {
340   assert(atom_u({ SuccEq({ ['A', true ], ['B', false] },
341                          { ['A', false], ['B', false] }) }) == { 'A', 'B' },
342          "Error: atom_u is incorrect.");
343   print('[ok] atom_u');
344 };
345 test_atom_u();
346 
347 test_card_u_to_ord_u := procedure() {
348   // In actual usage, instead of letters, states/events would be compared.
349   assert(card_u_to_ord_u({ ['A', 5], ['B', 3], ['C', 3], ['D', 1] },
350                          { 'A', 'B', 'C', 'D' }) ==
351          { SuccEq('A', 'B'), SuccEq('A', 'C'), SuccEq('A', 'D'),
352            SuccEq('B', 'C'), SuccEq('B', 'D'),
353            SuccEq('C', 'B'), SuccEq('C', 'D'),
354            // reflexive
355            SuccEq('A', 'A'), SuccEq('B', 'B'),
356            SuccEq('C', 'C'), SuccEq('D', 'D')
357          },
358          "Error: card_u_to_ord_u is incorrect.");
359   print('[ok] card_u_to_ord_u');
360 };
361 test_card_u_to_ord_u();
362 
363 test_add_u_to_card_u := procedure() {
364   assert(add_u_to_card_u({ ['A', 3], ['B', 5], ['C', 9],
365                            [And('A','B'), -2] }) ==
366          { [{ ['A', false], ['B', false], ['C', false] }, 0],
367            [{ ['A', false], ['B', false], ['C', true ] }, 9],
368            [{ ['A', false], ['B', true ], ['C', false] }, 5],
369            [{ ['A', false], ['B', true ], ['C', true ] }, 14],
370            [{ ['A', true ], ['B', false], ['C', false] }, 3],
371            [{ ['A', true ], ['B', false], ['C', true ] }, 12],
372            [{ ['A', true ], ['B', true ], ['C', false] }, 6],
373            [{ ['A', true ], ['B', true ], ['C', true ] }, 15]
374          },
375          "Error: add_u_to_card_u is incorrect.");
376   print('[ok] add_u_to_card_u');
377 };
378 test_add_u_to_card_u();
379 
380 test_ord_u_dist := procedure() {
381   dom := exprs_to_state_space({ 'A', 'B' });
382   ldom := [ x : x in dom ];
383   u1 := card_u_to_ord_u({
384     [ldom[1], 10], [ldom[2], 7], [ldom[3], 3 ], [ldom[4], 0]
385   }, dom);
386   u2 := card_u_to_ord_u({
387     [ldom[1], 9], [ldom[2], 1], [ldom[3], 11 ], [ldom[4], 4]
388   }, dom);
389   assert(ord_u_dist(u1, u2) == 1/2,
390          "Error: ord_u_dist is incorrect.");
391   print('[ok] ord_u_dist');
392   print('[  ] ord_u_dist w/ sharpenings');
393 };
394 test_ord_u_dist();
395 
396 test_sort_list := procedure() {
397   less := procedure(x, y) { return x < y; };
398   greater := procedure(x, y) { return y < x; };
399 
400   l := [1,3,5,4,2];
401   s1 := sort_list(l, less );
402   s2 := sort_list(l, greater);
403   assert(s1 == [1,2,3,4,5], 'Error: sort_list is incorrect.');
404   assert(s2 == [5,4,3,2,1], 'Error: sort_list is incorrect.');
405   print('[ok] sort_list');
406 };
407 test_sort_list();
408 
409 test_sort_list := procedure() {
410   less := procedure(x, y) { return x < y; };
411   s := {1,3,5,4,2};
412   assert(sort_set(s, less) == [1,2,3,4,5], 'Error: sort_list is incorrect.');
413   print('[ok] sort_set');
414 };
415 test_sort_list();
416 
417 test_restrict_domain := procedure() {
418   assert(restrict_domain({ ['a', 1], ['b', 3], ['c', 5] }, {'b', 'c'})
419            == { ['b', 3], ['c', 5] },
420          "Error: is incorrect.");
421   print('[ok] restrict_doman');
422 };
423 test_restrict_domain();
424 
425 test_restricted_eq := procedure() {
426   assert(restricted_eq({['a', 1], ['b', 2], ['c', 3]},
427                                  {['b', 2], ['c', 3]} ) == true,
428          "Error: restricted_eq true is incorrect.");
429   print('[ok] restricted_eq true');
430   assert(restricted_eq({['a', 1], ['b', 2], ['c', 3]},
431                                  {['b', 2], ['c', 3], ['d', 4]} ) == false,
432          "Error: restricted_eq false is incorrect.");
433   print('[ok] restricted_eq false');
434 };
435 test_restricted_eq();
436 
437 test_range_for_arity := procedure() {
438   poss_sts := possible_states({ 'A', 'B' }, { ['A', tf()], ['B', tf()] });
439   sts := { first(poss_sts), { ['A', true], ['B', false]} };
440   assert(sts in range_for_arity(0, poss_sts),
441          "Error: range_for_arity 0 is incorrect.");
442   print('[ok] range_for_arity 0');
443   sts2 := { { ['A', false], ['B', false]}, { ['A', false], ['B', true]} };
444   // assert({ [sts], [sts2] } in range_for_arity(1, poss_sts), 
445   //       "Error: range_for_arity 1 is incorrect.");
446   // print('[ok] range_for_arity 1');
447   print('[..] range_for_arity 1');
448   // assert({ [sts, sts2] } in range_for_arity(2, poss_sts),
449   //       "Error: range_for_arity 2 is incorrect.");
450   // print('[ok] range_for_arity 2');
451 };
452 test_range_for_arity();
453 
454 test_exprs_to_state_space := procedure() {
455   assert(exprs_to_state_space({ "A", "B", Implies("A", "B"),
456                                 Not(And("A", Not("B"))) }) == {
457            {['A', false], ['B', false]}, {['A', false], ['B', true]},
458            {['A',  true], ['B', false]}, {['A',  true], ['B', true]}
459          }, "Error: exprs_to_state_space is incorrect.");
460   print('[ok] exprs_to_state_space');
461 };
462 test_exprs_to_state_space();
1 /* SetlX 2.2.0 Cheat Sheet 2 Here I've compiled a short list of examples illustrating the setlX language. 3 4 There is also a slightly expanded quick_ref.stlx version which demonstrates 5 nearly all the setlX needed in order to understand MetaEthical.AI. The main 6 difference is that this cheat sheet focuses on aspects that are more unique 7 to the setlX language and omits many commands that are more obvious or 8 self-explanatory, especially among programmers. 9 10 A much more thorough guide can be found in the official setlX tutorial.pdf. 11 */ 12 13 // Standard Arithmetic Operations 14 x := (1 + 2) * 3 / 4; 15 // ~< Result: 9/4 >~ 16 /* Note that lines always end with a semi-colon. 17 Assignment of a value to a variable is done with ":=", not "=". 18 Variable names start with a lowercase letter because capitalized names are 19 reserved for setlx functors. */ 20 21 // Exponents 22 2 ** 5; 23 // ~< Result: 32 >~ 24 25 // String Length 26 #"abc"; 27 // ~< Result: 3 >~ 28 29 30 // Sets 31 // ==== 32 a := { 1, 2, 3 }; 33 b := { 3, 4, 5 }; 34 35 // Cardinality 36 #b; 37 // ~< Result: 3 >~ 38 39 // Union 40 a + b; 41 // ~< Result: {1, 2, 3, 4, 5} >~ 42 43 // Intersection 44 a * b; 45 // ~< Result: {3} >~ 46 47 // Difference 48 a - b; 49 // ~< Result: {1, 2} >~ 50 51 // Check Membership 52 1 in a; 53 // ~< Result: true >~ 54 55 // Check if Subset 56 {1, 2} < {1, 2, 3}; 57 // ~< Result: true >~ 58 59 // Powerset 60 pow({ 1, 2 }); 61 // ~< Result: {{}, {1}, {1, 2}, {2}} >~ 62 2 ** { 1, 2 }; // alternate syntax 63 // ~< Result: {{}, {1}, {1, 2}, {2}} >~ 64 65 // Cartesian Product 66 pairs := { 'a', 'b' } >< { 1, 2 }; 67 // ~< Result: {["a", 1], ["a", 2], ["b", 1], ["b", 2]} >~ 68 69 // Set comprehension 70 { x * 2 : x in a }; 71 // ~< Result: {2, 4, 6} >~ 72 73 // Set comprehension with filter 74 { x : x in a | x < 3 }; 75 // ~< Result: {1, 2} >~ 76 77 // Set comprehension with list pattern matching 78 // Also, since we don't use the first element of the ordered pair, we can 79 // discard it with _ instead of giving it a variable name. 80 { y : [_, y] in pairs }; 81 // ~< Result: {1, 2} >~ 82 83 // Quantification 84 exists(x in a | x == 2); 85 // ~< Result: true >~ 86 forall(x in a | x < 5); 87 // ~< Result: true >~ 88 89 // Returns the first element according to the language's internal order. 90 // If the set is not all numbers or lists of them, the ordering should not be 91 // relied upon. 92 first(a); 93 // ~< Result: 1 >~ 94 95 // Sum of a Set of... 96 // Integers 97 +/ {1, 2, 3}; 98 // ~< Result: 6 >~ 99 100 // Strings 101 +/ {'a', 'b', 'c'}; 102 // ~< Result: "abc" >~ 103 104 // Sets 105 +/ { {1, 2}, {3, 4} }; 106 // ~< Result: {1, 2, 3, 4} >~ 107 108 // Lists 109 +/ { [1, 2], [3, 4] }; 110 // ~< Result: [1, 2, 3, 4] >~ 111 112 113 // Lists 114 // ===== 115 l := [ 1, 2, 3 ]; 116 m := [ 3, 4, 5 ]; 117 118 // Note that list indexes start at 1! 119 // Not 0 as with most other programming languages. 120 m[1]; 121 // ~< Result: 3 >~ 122 123 // Concatenation 124 l + m; 125 // ~< Result: [ 1, 2, 3, 3, 4, 5] >~ 126 // Note that unlike with sets, the 3 is duplicated 127 128 // List comprehension 129 [ x * 2 : x in l ]; 130 // ~< Result: [2, 4, 6] >~ 131 132 // List comprehension with filter 133 [ x : x in l | x < 3 ]; 134 // ~< Result: [1, 2] >~ 135 136 // Check Inclusion 137 1 in l; 138 // ~< Result: true >~ 139 140 // Quantification 141 exists(x in l | x == 2); 142 // ~< Result: true >~ 143 forall(x in l | x < 5); 144 // ~< Result: true >~ 145 146 // Range 147 m[1..2]; 148 // ~< Result: [3, 4] >~ 149 150 m[..2]; 151 // ~< Result: [3, 4] >~ 152 153 m[2..]; 154 // ~< Result: [4, 5] >~ 155 156 first([3,2,1]); 157 // ~< Result: 3 >~ 158 159 last([3,2,1]); 160 // ~< Result: 1 >~ 161 162 permutations([1, 2, 3]); 163 // ~< Result: {[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]} >~ 164 165 // Sum of a List of: 166 // Integers 167 +/ [1, 2, 3]; 168 // ~< Result: 6 >~ 169 170 // Strings 171 +/ ['a', 'b', 'c']; 172 // ~< Result: "abc" >~ 173 174 // Sets 175 +/ [ {1, 2}, {3, 4} ]; 176 // ~< Result: {1, 2, 3, 4} >~ 177 178 // Lists 179 +/ [ [1, 2], [3, 4] ]; 180 // ~< Result: [1, 2, 3, 4] >~ 181 182 183 // Functions 184 // ========= 185 // They are just sets of ordered pairs where each element of the domain is 186 // mapped to one element of the range. 187 f := { [1,2], [2,4], [3,6] }; 188 189 // You can find the value given an input 190 f[2]; 191 // ~< Result: 4 >~ 192 193 // But if it's not a function (an input is mapped to more than one output) 194 // then it's om (which is like javascript's null). 195 g := { [1,2], [1,3] }; 196 g[1]; 197 // ~< Result: om >~ 198 199 domain(f); 200 // ~< Result: {1, 2, 3} >~ 201 202 range(f); 203 // ~< Result: {2, 4, 6} >~ 204 205 isMap({ [1,2], [2,4] }); 206 // ~< Result: true >~ 207 isMap({ [1,2], [1,4] }); 208 // ~< Result: false >~ 209 210 211 // Control Structures 212 // ================== 213 switch { 214 case false : 215 print('hi'); 216 case false : 217 print('hey'); 218 default : 219 print('bye'); 220 } 221 // bye 222 223 224 // Procedures 225 sum := procedure(x, y) { 226 return x + y; 227 }; 228 sum(3, 4); 229 // ~< Result: 7 >~ 230 231 // It can also be cached to avoid recomputing when given the same inputs 232 sum := cachedProcedure(x, y) { 233 return x + y; 234 }; 235 sum(3, 4); 236 // ~< Result: 7 >~ 237 238 239 // Object-oriented Programming 240 class point(x, y) { 241 x := x; 242 y := y; 243 244 // Instance procedure 245 dist_to := procedure(pt) { 246 return sqrt( (pt.x - x) ** 2 + (pt.y - get_y()) ** 2 ); 247 }; 248 249 // Not really necessary but just demonstrating calling another instance proc 250 get_y := procedure() { 251 // Also just demonstrating 'this' keyword, usually unneeded and implicit 252 return this.y; 253 }; 254 255 // Class procedures 256 static { 257 orig := procedure() { 258 return point(0, 0); 259 }; 260 } 261 } 262 pt := point(1, 1); 263 pt2 := point(4, 5); 264 pt.dist_to(pt2); 265 // ~< Result: 5.0 >~ 266 pt3 := point(3, 4); 267 pt3.dist_to(point.orig()); 268 // ~< Result: 5.0 >~ 269 270 271 /* SetlX Functors 272 ============== 273 SetlX allows for the creation of data structures composed from what they call 274 functors and arguments. Functor names always start with a capital letter and 275 then are followed by parentheses enclosing its arguments, which can be of 276 more or less any type. And that's basically all they are, like strings that 277 can take arguments to enable tree-like structures. 278 279 Here we demonstrate its usage to represent logical expressions where 280 arguments are either strings or other functor terms. (If you're familiar 281 with functors from category theory, the usage of that term here doesn't 282 really have anything to do with that notion.) 283 */ 284 And('A', 'B'); 285 // ~< Result: And("A", "B") >~ 286 287 Implies(And('A', 'B'), 'C'); 288 // ~< Result: Implies(And("A", "B"), "C") >~ 289 290 fct(And('A', 'B')); 291 // ~< Result: "And" >~ 292 293 args(And('A', 'B')); 294 // ~< Result: ["A", "B"] >~ 295 296 makeTerm('And', ['A', 'B']); 297 // ~< Result: And("A", "B") >~ 298 299 300 // Meta-programming 301 eval("{1} + {2, 3}"); 302 // ~< Result: {1, 2, 3} >~
1 /* SetlX 2.2.0 Quick Reference 2 This list of examples introduces much of the setlX language. 3 4 Nearly all of MetaEthical.AI is just a combination of the basic components 5 demonstrated here. There is also a slightly more compact cheat_sheet.stlx 6 that focuses on just the aspects that are more unique to setlX and omits 7 many commands that are more obvious or self-explanatory, especially among 8 programmers. 9 10 A much more thorough guide can be found in the official setlX tutorial.pdf. 11 */ 12 13 14 // This is a comment. 15 /* This is a 16 multi-line comment */ 17 18 19 // Standard Arithmetic Operations 20 x := (1 + 2) * 3 / 4; 21 // ~< Result: 9/4 >~ 22 /* Note that lines always end with a semi-colon. 23 Assignment of a value to a variable is done with ":=", not "=". 24 Variable names start with a lowercase letter because capitalized names are 25 reserved for setlx functors. */ 26 27 // Exponents 28 2 ** 5; 29 // ~< Result: 32 >~ 30 31 sqrt(81); 32 // ~< Result: 9.0 >~ 33 34 // Absolute value 35 abs(-7); 36 // ~< Result: 7 >~ 37 38 39 // Strings 40 // ======= 41 str1 := "abc"; 42 str2 := 'def'; // single quotes work too 43 44 // String concatenation 45 str3 := str1 + str2; 46 // ~< Result: "abcdef" >~ 47 48 // Length 49 #str3; 50 // ~< Result: 6 >~ 51 52 53 // Boolean values 54 // ============== 55 p := true; 56 q := false; 57 58 // Disjunction 59 p || q; 60 // ~< Result: true >~ 61 62 // Conjunction 63 p && q; 64 // ~< Result: false >~ 65 66 // Equality / Biconditional 67 p == q; 68 // ~< Result: false >~ 69 70 p != q; 71 // ~< Result: true >~ 72 73 // Negation 74 !p; 75 // ~< Result: false >~ 76 77 78 // Sets 79 // ==== 80 a := { 1, 2, 3 }; 81 b := { 3, 4, 5 }; 82 83 // Cardinality 84 #b; 85 // ~< Result: 3 >~ 86 87 // Union 88 a + b; 89 // ~< Result: {1, 2, 3, 4, 5} >~ 90 91 // Intersection 92 a * b; 93 // ~< Result: {3} >~ 94 95 // Difference 96 a - b; 97 // ~< Result: {1, 2} >~ 98 99 // Check Membership 100 1 in a; 101 // ~< Result: true >~ 102 103 // Check if Subset 104 {1, 2} < {1, 2, 3}; 105 // ~< Result: true >~ 106 {3} < {1, 2}; 107 // ~< Result: false >~ 108 109 // Powerset 110 pow({ 1, 2 }); 111 // ~< Result: {{}, {1}, {1, 2}, {2}} >~ 112 2 ** { 1, 2 }; // alternate syntax 113 // ~< Result: {{}, {1}, {1, 2}, {2}} >~ 114 115 // Cartesian Product 116 pairs := { 'a', 'b' } >< { 1, 2 }; 117 // ~< Result: {["a", 1], ["a", 2], ["b", 1], ["b", 2]} >~ 118 119 // Set comprehension 120 { x * 2 : x in a }; 121 // ~< Result: {2, 4, 6} >~ 122 123 // Set comprehension with filter 124 { x : x in a | x < 3 }; 125 // ~< Result: {1, 2} >~ 126 127 // Set comprehension with list pattern matching 128 // Also, since we don't use the first element of the ordered pair, we can 129 // discard it with _ instead of giving it a variable name. 130 { y : [_, y] in pairs }; 131 // ~< Result: {1, 2} >~ 132 133 // Quantification 134 exists(x in a | x == 2); 135 // ~< Result: true >~ 136 forall(x in a | x < 5); 137 // ~< Result: true >~ 138 139 // Returns the first element according to the language's internal order. 140 // If the set is not all numbers or lists of them, the ordering should not be 141 // relied upon. 142 first(a); 143 // ~< Result: 1 >~ 144 145 // Sum of a Set of... 146 // Integers 147 +/ {1, 2, 3}; 148 // ~< Result: 6 >~ 149 150 // Strings 151 +/ {'a', 'b', 'c'}; 152 // ~< Result: "abc" >~ 153 154 // Sets 155 +/ { {1, 2}, {3, 4} }; 156 // ~< Result: {1, 2, 3, 4} >~ 157 158 // Lists 159 +/ { [1, 2], [3, 4] }; 160 // ~< Result: [1, 2, 3, 4] >~ 161 162 163 // Lists 164 // ===== 165 l := [ 1, 2, 3 ]; 166 m := [ 3, 4, 5 ]; 167 168 // Note that list indexes start at 1! 169 // Not 0 as with most other programming languages. 170 m[1]; 171 // ~< Result: 3 >~ 172 173 // Concatenation 174 l + m; 175 // ~< Result: [ 1, 2, 3, 3, 4, 5] >~ 176 // Note that unlike with sets, the 3 is duplicated 177 178 // List comprehension 179 [ x * 2 : x in l ]; 180 // ~< Result: [2, 4, 6] >~ 181 182 // List comprehension with filter 183 [ x : x in l | x < 3 ]; 184 // ~< Result: [1, 2] >~ 185 186 // Check Inclusion 187 1 in l; 188 // ~< Result: true >~ 189 190 // Quantification 191 exists(x in l | x == 2); 192 // ~< Result: true >~ 193 forall(x in l | x < 5); 194 // ~< Result: true >~ 195 196 // Range 197 m[1..2]; 198 // ~< Result: [3, 4] >~ 199 200 m[..2]; 201 // ~< Result: [3, 4] >~ 202 203 m[2..]; 204 // ~< Result: [4, 5] >~ 205 206 first([3,2,1]); 207 // ~< Result: 3 >~ 208 209 last([3,2,1]); 210 // ~< Result: 1 >~ 211 212 permutations([1, 2, 3]); 213 // ~< Result: {[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]} >~ 214 215 // Sum of a List of: 216 // Integers 217 +/ [1, 2, 3]; 218 // ~< Result: 6 >~ 219 220 // Strings 221 +/ ['a', 'b', 'c']; 222 // ~< Result: "abc" >~ 223 224 // Sets 225 +/ [ {1, 2}, {3, 4} ]; 226 // ~< Result: {1, 2, 3, 4} >~ 227 228 // Lists 229 +/ [ [1, 2], [3, 4] ]; 230 // ~< Result: [1, 2, 3, 4] >~ 231 232 233 // Functions 234 // ========= 235 // They are just sets of ordered pairs where each element of the domain is 236 // mapped to one element of the range. 237 f := { [1,2], [2,4], [3,6] }; 238 239 // You can find the value given an input 240 f[2]; 241 // ~< Result: 4 >~ 242 243 // But if it's not a function (an input is mapped to more than one output) 244 // then it's om (which is like javascript's null). 245 g := { [1,2], [1,3] }; 246 g[1]; 247 // ~< Result: om >~ 248 249 domain(f); 250 // ~< Result: {1, 2, 3} >~ 251 252 range(f); 253 // ~< Result: {2, 4, 6} >~ 254 255 256 // Control Structures 257 // ================== 258 if (true) { 259 print('hi'); 260 } 261 // hi 262 263 if (false) { 264 print('hi'); 265 } else { 266 print('bye'); 267 } 268 // bye 269 270 if (false) { 271 print('hi'); 272 } else if (true) { 273 print('hey'); 274 } else { 275 print('bye'); 276 } 277 // hey 278 279 switch { 280 case false : 281 print('hi'); 282 case false : 283 print('hey'); 284 default : 285 print('bye'); 286 } 287 // bye 288 289 s := {}; 290 for (i in [1, 2, 3]) { 291 s += { i*2 }; 292 } 293 s; 294 // ~< Result: {2, 4, 6} >~ 295 296 x := 0; 297 while (x < 5) { 298 x += 1; 299 } 300 x; 301 // ~< Result: 5 >~ 302 303 // Todo: match, _ 304 305 306 // Procedures 307 sum := procedure(x, y) { 308 return x + y; 309 }; 310 sum(3, 4); 311 // ~< Result: 7 >~ 312 313 // It can also be cached to avoid recomputing when given the same inputs 314 sum := cachedProcedure(x, y) { 315 return x + y; 316 }; 317 sum(3, 4); 318 // ~< Result: 7 >~ 319 320 321 // Object-oriented Programming 322 class point(x, y) { 323 x := x; 324 y := y; 325 326 // Instance procedure 327 dist_to := procedure(pt) { 328 return sqrt( (pt.x - x) ** 2 + (pt.y - get_y()) ** 2 ); 329 }; 330 331 // Not really necessary but just demonstrating calling another instance proc 332 get_y := procedure() { 333 // Also just demonstrating 'this' keyword, usually unneeded and implicit 334 return this.y; 335 }; 336 337 // Class procedures 338 static { 339 orig := procedure() { 340 return point(0, 0); 341 }; 342 } 343 } 344 pt := point(1, 1); 345 pt2 := point(4, 5); 346 pt.dist_to(pt2); 347 // ~< Result: 5.0 >~ 348 pt3 := point(3, 4); 349 pt3.dist_to(point.orig()); 350 // ~< Result: 5.0 >~ 351 352 353 /* SetlX Functors 354 ============== 355 SetlX allows for the creation of data structures composed from what they call 356 functors and arguments. Functor names always start with a capital letter and 357 then are followed by parentheses enclosing its arguments, which can be of 358 more or less any type. And that's basically all they are, like strings that 359 can take arguments to enable tree-like structures. 360 361 Here we demonstrate its usage to represent logical expressions where 362 arguments are either strings or other functor terms. (If you're familiar 363 with functors from category theory, the usage of that term here doesn't 364 really have anything to do with that notion.) 365 */ 366 And('A', 'B'); 367 // ~< Result: And("A", "B") >~ 368 369 Implies(And('A', 'B'), 'C'); 370 // ~< Result: Implies(And("A", "B"), "C") >~ 371 372 fct(And('A', 'B')); 373 // ~< Result: "And" >~ 374 375 args(And('A', 'B')); 376 // ~< Result: ["A", "B"] >~ 377 378 makeTerm('And', ['A', 'B']); 379 // ~< Result: And("A", "B") >~ 380 381 382 // Type checking 383 // ============= 384 isTerm(And('A', 'B')); 385 // ~< Result: true >~ 386 387 isString('abc'); 388 // ~< Result: true >~ 389 390 isSet({}); 391 // ~< Result: true >~ 392 393 isList([]); 394 // ~< Result: true >~ 395 396 isMap({ [1,2], [2,4] }); 397 // ~< Result: true >~ 398 isMap({ [1,2], [1,4] }); 399 // ~< Result: false >~ 400 401 402 int('5'); 403 // ~< Result: 5 >~ 404 405 str(5); 406 // ~< Result: "5" >~ 407 408 str({}); 409 // ~< Result: "{}" >~ 410 411 412 // Meta-programming 413 eval("{1} + {2, 3}"); 414 // ~< Result: {1, 2, 3} >~