*

3 files changed, 649 insertions, 0 deletions

diff --git a/picat/planner.pi b/picat/planner.pi
new file mode 100644
index 0000000..8cca7b2
--- /dev/null
+++ b/picat/planner.pi
@@ -0,0 +1,320 @@
+ module planner.
+
+% best_plan(S,Limit,Plan) => best_plan(S,Limit,Plan).
+% best_plan(S,Limit,Plan,PlanCost) => best_plan(S,Limit,Plan,PlanCost).
+% best_plan(S,Plan,PlanCost) => best_plan(S,Plan,PlanCost).
+% best_plan(S,Plan) => best_plan(S,Plan).
+% best_plan_bin(S,Limit,Plan) => best_plan_bin(S,Limit,Plan).
+% best_plan_bin(S,Limit,Plan,PlanCost) => best_plan_bin(S,Limit,Plan,PlanCost).
+% best_plan_bin(S,Plan,PlanCost) => best_plan_bin(S,Plan,PlanCost).
+% best_plan_bin(S,Plan) => best_plan_bin(S,Plan).
+% best_plan_bb(S,Limit,Plan) => best_plan_bb(S,Limit,Plan).
+% best_plan_bb(S,Plan,PlanCost) => best_plan_bb(S,Plan,PlanCost).
+% best_plan_bb(S,Limit,Plan,PlanCost) => best_plan_bb(S,Limit,Plan,PlanCost).
+% best_plan_bb(S,Plan) => best_plan_bb(S,Plan).
+% best_plan_nondet(S,Limit,Plan) => best_plan_nondet(S,Limit,Plan).
+% best_plan_nondet(S,Plan,PlanCost) => best_plan_nondet(S,Plan,PlanCost).
+% best_plan_nondet(S,Limit,Plan,PlanCost) => best_plan_nondet(S,Limit,Plan,PlanCost).
+% best_plan_nondet(S,Plan) => best_plan_nondet(S,Plan).
+% best_plan_unbounded(S,Limit,Plan) => best_plan_unbounded(S,Limit,Plan).
+% best_plan_unbounded(S,Plan,PlanCost) => best_plan_unbounded(S,Plan,PlanCost).
+% best_plan_unbounded(S,Limit,Plan,PlanCost) => best_plan_unbounded(S,Limit,Plan,PlanCost).
+% best_plan_unbounded(S,Plan) => best_plan_unbounded(S,Plan).
+% current_plan() = current_plan().
+% current_resource() = current_resource().
+% current_resource_plan(Amount,Plan) => current_resource_plan(Amount,Plan).
+% current_resource_plan_cost(Amount,Plan,Cost) => current_resource_plan_cost(Amount,Plan,Cost).
+% insert_state_list(StateL,Elm) = insert_state_list(StateL,Elm).
+% is_tabled_state(S) => is_tabled_state(S).
+% new_state_list(List) = new_state_list(List).
+% plan(S,Limit,Plan) => plan(S,Limit,Plan).
+% plan(S,Plan,PlanCost) => plan(S,Plan,PlanCost).
+% plan(S,Limit,Plan,PlanCost) => plan(S,Limit,Plan,PlanCost).
+% plan(S,Plan) => plan(S,Plan).
+% plan_unbounded(S,Limit,Plan) => plan_unbounded(S,Limit,Plan).
+% plan_unbounded(S,Plan,PlanCost) => plan_unbounded(S,Plan,PlanCost).
+% plan_unbounded(S,Limit,Plan,PlanCost) => plan_unbounded(S,Limit,Plan,PlanCost).
+% plan_unbounded(S,Plan) => plan_unbounded(S,Plan).
+
+% A state-list is an ordered list, where the ordering of symbols is determined 
+% by the addresses of the symbols in the symbol table, not by the characters that
+% constitute the symbols. This ordering allows for efficient ordering of symbols. 
+% Note that this ordering is different from lexicographical ordering, which is used by sort(). 
+
+new_state_list(List) = SList => new_state_list_aux(List,[],SList).
+
+new_state_list_aux([],SList0,SList) => SList = SList0.
+new_state_list_aux([E|List],SList0,SList) => 
+    bp.b_INSERT_STATE_LIST_ccf(SList0,E,SList1),
+    new_state_list_aux(List,SList1,SList).
+
+insert_state_list(SList0, Elm) = SList => bp.b_INSERT_STATE_LIST_ccf(SList0,Elm,SList).
+
+%%%
+current_resource() = Amount =>
+    bp.b_PLANNER_CURR_RPC_fff(Amount,_Plan,_Cost).
+
+current_plan() = Plan =>
+    bp.b_PLANNER_CURR_RPC_fff(_Amount,Plan,_Cost).
+
+current_resource_plan(Amount,Plan) =>
+    bp.b_PLANNER_CURR_RPC_fff(Amount,Plan,_Cost).
+
+current_resource_plan_cost(Amount,Plan,Cost) =>
+    bp.b_PLANNER_CURR_RPC_fff(Amount,Plan,Cost).
+
+%%%
+plan(S,Plan), var(Plan) =>
+    bp.picat_plan(S,268435455,Plan,_).    % the limit is assumed to be 268435455
+plan(_S,Plan) =>
+    throw_plan_arg_error(1,Plan,_,plan).
+
+plan(S,Limit,Plan), var(Plan), integer(Limit), Limit >= 0 =>
+    bp.picat_plan(S,Limit,Plan,_). 
+plan(S,Plan,PlanCost), var(Plan), var(PlanCost) =>
+    bp.picat_plan(S,268435455,Plan,PlanCost). 
+plan(_S,Limit,Plan) =>
+    throw_plan_arg_error(Limit,Plan,_,plan).
+
+plan(S,Limit,Plan,PlanCost), var(Plan), var(PlanCost), integer(Limit), Limit >= 0 =>
+    bp.picat_plan(S,Limit,Plan,PlanCost). 
+plan(_S,Limit,Plan,PlanCost) =>
+    throw_plan_arg_error(Limit,Plan,PlanCost,plan).
+
+plan(S,Limit,Plan,PlanCost,FinS), var(Plan), var(PlanCost), integer(Limit), Limit >= 0 =>
+    bp.picat_plan(S,Limit,Plan,PlanCost,FinS). 
+plan(_S,Limit,Plan,PlanCost,_FinS) =>
+    throw_plan_arg_error(Limit,Plan,PlanCost,plan).
+
+%%%
+plan_unbounded(S,Plan), var(Plan) =>
+    bp.picat_plan_unbounded(S,268435455,Plan,_).    % the limit is assumed to be 268435455
+plan_unbounded(_S,Plan) =>
+    throw_plan_arg_error(1,Plan,_,plan_unbounded).
+
+plan_unbounded(S,Limit,Plan), var(Plan), integer(Limit), Limit >= 0 =>
+    bp.picat_plan_unbounded(S,Limit,Plan,_). 
+plan_unbounded(S,Plan,PlanCost), var(Plan), var(PlanCost) =>
+    bp.picat_plan_unbounded(S,268435455,Plan,PlanCost). 
+plan_unbounded(_S,Limit,Plan) =>
+    throw_plan_arg_error(Limit,Plan,_,plan_unbounded).
+
+plan_unbounded(S,Limit,Plan,PlanCost), var(Plan), var(PlanCost), integer(Limit), Limit >= 0 =>
+    bp.picat_plan_unbounded(S,Limit,Plan,PlanCost). 
+plan_unbounded(_S,Limit,Plan,PlanCost) =>
+    throw_plan_arg_error(Limit,Plan,PlanCost,plan_unbounded).
+
+
+%%%
+%%% iterative deepening
+best_plan(S,Plan), var(Plan) =>
+    best_plan_downward(S,0,268435455,Plan,_,_).
+best_plan(_S,Plan) =>
+    throw_plan_arg_error(1,Plan,_,best_plan).
+
+best_plan(S,Limit,Plan), var(Plan), integer(Limit), Limit >= 0 =>
+    best_plan_downward(S,0,Limit,Plan,_,_).
+best_plan(S,Plan,PlanCost), var(Plan), var(PlanCost) =>
+    best_plan_downward(S,0,268435455,Plan,PlanCost,_).
+best_plan(_S,Limit,Plan) =>
+    throw_plan_arg_error(Limit,Plan,_,best_plan).
+
+best_plan(S,Limit,Plan,PlanCost), var(Plan), var(PlanCost), integer(Limit), Limit >= 0 =>
+    best_plan_downward(S,0,Limit,Plan,PlanCost,_).
+best_plan(_S,Limit,Plan,PlanCost) =>
+    throw_plan_arg_error(Limit,Plan,PlanCost,best_plan).
+
+best_plan(S,Limit,Plan,PlanCost,FinS), var(Plan), var(PlanCost), integer(Limit), Limit >= 0 =>
+    best_plan_downward(S,0,Limit,Plan,PlanCost,FinS).
+best_plan(_S,Limit,Plan,PlanCost,_FinS) =>
+    throw_plan_arg_error(Limit,Plan,PlanCost,best_plan).
+
+%%%
+best_plan_downward(S,Level,_Limit,_Plan,_PlanCost,_FinS) ?=>
+    (bp.global_get('_$picat_log',0, 1) -> printf("%% Searching with the bound %d\n",Level); true),
+    call_picat_plan(S,Level),
+    fail.
+best_plan_downward(S,_Level,_Limit,Plan,PlanCost,FinS),
+    Map = get_table_map('_$planner'),
+    Map.has_key($current_best_plan(S))
+ =>
+    Map.get($current_best_plan(S)) = Plan,
+    Map.get($current_best_plan_cost(S)) = PlanCost,
+    Map.get($current_best_plan_fin_state(S)) = FinS.
+best_plan_downward(S,Level,Limit,Plan,PlanCost,FinS) =>
+    bp.global_get('_$planner_explored_depth',0,Depth),
+    (Depth == 268435455 -> Level1 = Level+1; Level1 = Depth),
+    Level1 =< Limit,
+    best_plan_downward(S,Level1,Limit,Plan,PlanCost,FinS).
+
+%%% iterative deeping and binary search
+best_plan_bin(S,Plan), var(Plan) =>
+    best_plan_downward_bin(S,0,268435455,Plan,_,_).
+best_plan_bin(_S,Plan) =>
+    throw_plan_arg_error(1,Plan,_,best_plan).
+
+best_plan_bin(S,Limit,Plan), var(Plan), integer(Limit), Limit >= 0 =>
+    best_plan_downward_bin(S,0,Limit,Plan,_,_).
+best_plan_bin(S,Plan,PlanCost), var(Plan), var(PlanCost) =>
+    best_plan_downward_bin(S,0,268435455,Plan,PlanCost,_).
+best_plan_bin(_S,Limit,Plan) =>
+    throw_plan_arg_error(Limit,Plan,_,best_plan).
+
+best_plan_bin(S,Limit,Plan,PlanCost), var(Plan), var(PlanCost), integer(Limit), Limit >= 0 =>
+    best_plan_downward_bin(S,0,Limit,Plan,PlanCost,_).
+best_plan_bin(_S,Limit,Plan,PlanCost) =>
+    throw_plan_arg_error(Limit,Plan,PlanCost,best_plan).
+
+best_plan_bin(S,Limit,Plan,PlanCost,FinS), var(Plan), var(PlanCost), integer(Limit), Limit >= 0 =>
+    best_plan_downward_bin(S,0,Limit,Plan,PlanCost,FinS).
+best_plan_bin(_S,Limit,Plan,PlanCost,_FinS) =>
+    throw_plan_arg_error(Limit,Plan,PlanCost,best_plan).
+
+%%%
+best_plan_downward_bin(S,Level,Limit,Plan,PlanCost,FinS) =>
+    Map = get_table_map('_$planner'),
+    loop_best_plan_downward_bin(S,Level,Limit,Plan,PlanCost,FinS,Map).
+
+loop_best_plan_downward_bin(S,Level,Limit,_Plan,_PlanCost,_FinS,_Map) ?=>
+    Level =< Limit,
+    (bp.global_get('_$picat_log',0, 1) -> printf("%% Searching with the bound %d\n",Level); true),
+    call_picat_plan(S,Level),
+    fail.
+loop_best_plan_downward_bin(S,_Level,_Limit,Plan,PlanCost,_FinS,Map),
+    Map.has_key($current_best_plan(S))
+ =>
+    Lower = Map.get($current_lower_bound(S), 0),
+    Upper = Map.get($current_best_plan_cost(S))-1,
+    loop_best_plan_bin(S,Map,Lower,Upper,Plan,PlanCost).
+loop_best_plan_downward_bin(S,Level,Limit,Plan,PlanCost,FinS,Map) =>
+    bp.global_get('_$planner_explored_depth',0,Depth),
+%    writeln(depth=Depth),
+    (Depth == 268435455 -> Lower = Level+1; Lower = Depth),
+    Map.put($current_lower_bound(S), Lower),
+    Lower < Limit,
+    Level1 = 2*Lower,
+    NewLevel = cond(Level1 > Limit, Limit, Level1),
+    loop_best_plan_downward_bin(S,NewLevel,Limit,Plan,PlanCost,FinS,Map).
+
+% binary search
+loop_best_plan_bin(S,Map,Lower,Upper,Plan,PlanCost),
+    Lower =< Upper
+ =>
+    NewLimit = Lower + (Upper-Lower) div 2,
+    (bp.global_get('_$picat_log',0, 1) -> printf("%% Searching with the bound %d\n",NewLimit); true),
+    if call_picat_plan(S,NewLimit) then
+        NewUpper = Map.get($current_best_plan_cost(S))-1,
+        loop_best_plan_bin(S,Map,Lower,NewUpper,Plan,PlanCost)
+    else
+        NewLower = NewLimit+1,
+        loop_best_plan_bin(S,Map,NewLower,Upper,Plan,PlanCost)
+    end.
+loop_best_plan_bin(S,Map,_Lower,_Upper,Plan,PlanCost) =>
+    Map.has_key($current_best_plan(S)),
+    Map.get($current_best_plan(S)) = Plan,
+    Map.get($current_best_plan_cost(S)) = PlanCost.
+
+%%%
+best_plan_nondet(S,Plan), var(Plan) =>
+    best_plan_nondet_aux(S,268435455,Plan,_).
+best_plan_nondet(_S,Plan) =>
+    throw_plan_arg_error(1,Plan,_,best_plan_nondet).
+
+best_plan_nondet(S,Limit,Plan), var(Plan), integer(Limit), Limit >= 0 =>
+    best_plan_nondet_aux(S,Limit,Plan,_).
+best_plan_nondet(S,Plan,PlanCost), var(Plan), var(PlanCost) =>
+    best_plan_nondet_aux(S,268435455,Plan,PlanCost).
+best_plan_nondet(_S,Limit,Plan) =>
+    throw_plan_arg_error(Limit,Plan,_,best_plan_nondet).
+
+best_plan_nondet(S,Limit,Plan,PlanCost), var(Plan), var(PlanCost), integer(Limit), Limit >= 0 =>
+    best_plan_nondet_aux(S,Limit,Plan,PlanCost).
+best_plan_nondet(_S,Limit,Plan,PlanCost) =>
+    throw_plan_arg_error(Limit,Plan,PlanCost,best_plan_nondet).
+
+best_plan_nondet_aux(S,Limit,Plan,PlanCost) =>
+    not not (M = get_global_map(),
+             M.put('_first_best_plan',[]),
+             best_plan_downward(S,0,Limit,Plan0,PlanCost0,_),     % use tabled search to find the first best plan
+             M.put('_first_best_plan', (Plan0,PlanCost0))),
+    get_global_map().get('_first_best_plan') = (Plan0,PlanCost),
+    (    Plan = Plan0
+    ;
+        bp.picat_best_plan_nondet_nontabled_search(S,Plan,PlanCost),
+        Plan0 != Plan
+    ).
+
+%%% Branch-and-Bound
+best_plan_bb(S,Plan), var(Plan) =>
+    loop_best_plan_bb(S,268435455,Plan,_).
+best_plan_bb(_S,Plan) =>
+    throw_plan_arg_error(1,Plan,_,best_plan_bb).
+
+best_plan_bb(S,Limit,Plan), var(Plan), integer(Limit), Limit >= 0 =>
+    loop_best_plan_bb(S,Limit,Plan,_PlanCost).
+best_plan_bb(S,Plan,PlanCost), var(Plan), var(PlanCost) =>
+    loop_best_plan_bb(S,268435455,Plan,PlanCost).
+best_plan_bb(_S,Limit,Plan) =>
+    throw_plan_arg_error(Limit,Plan,_,best_plan_bb).
+
+best_plan_bb(S,Limit,Plan,PlanCost), var(Plan), var(PlanCost), integer(Limit), Limit >= 0 =>
+    loop_best_plan_bb(S,Limit,Plan,PlanCost).
+best_plan_bb(_S,Limit,Plan,PlanCost) =>
+    throw_plan_arg_error(Limit,Plan,PlanCost,best_plan_bb).
+
+loop_best_plan_bb(S,Limit,Plan,PlanCost),
+    (bp.global_get('_$picat_log',0, 1) -> printf("%% Searching with the bound %d\n",Limit); true),
+    call_picat_plan(S,Limit)        
+ =>
+    get_table_map('_$planner').get($current_best_plan_cost(S)) = PlanCost1,    
+    loop_best_plan_bb(S,PlanCost1-1,Plan,PlanCost).
+loop_best_plan_bb(S,_Limit,Plan,PlanCost) =>
+    Map = get_table_map('_$planner'),
+    Map.has_key($current_best_plan(S)),
+    Map.get($current_best_plan(S)) = Plan,
+    Map.get($current_best_plan_cost(S)) = PlanCost.
+
+call_picat_plan(S,Limit) =>
+    bp.global_set('_$planner_explored_depth',0,268435455),
+    not not call_picat_plan_aux(S,Limit).   % discard exception catchers created by picat_plan
+
+call_picat_plan_aux(S,Limit) =>
+    bp.picat_plan(S,Limit,Plan,PlanCost,FinS),
+    Map = get_table_map('_$planner'),
+    Map.put($current_best_plan(S), Plan),
+    Map.put($current_best_plan_cost(S), PlanCost),
+    Map.put($current_best_plan_fin_state(S), FinS).
+
+%%%
+best_plan_unbounded(S,Plan), var(Plan) =>
+    bp.picat_best_plan_unbounded(S,Plan,_).
+best_plan_unbounded(_S,Plan) =>
+    throw_plan_arg_error(1,Plan,_,best_plan_unbounded).
+
+best_plan_unbounded(S,Limit,Plan), var(Plan), integer(Limit), Limit >= 0 =>
+    bp.picat_best_plan_unbounded(S,Plan,PlanCost),
+    PlanCost =< Limit.
+best_plan_unbounded(S,Plan,PlanCost), var(Plan), var(PlanCost) =>
+    bp.picat_best_plan_unbounded(S,Plan,PlanCost).
+best_plan_unbounded(_S,Limit,Plan) =>
+    throw_plan_arg_error(Limit,Plan,_,best_plan_unbounded).    
+
+best_plan_unbounded(S,Limit,Plan,PlanCost), var(Plan), var(PlanCost), integer(Limit), Limit >= 0 =>
+    bp.picat_best_plan_unbounded(S,Plan,PlanCost),
+    PlanCost =< Limit.
+best_plan_unbounded(_S,Limit,Plan,PlanCost) =>
+    throw_plan_arg_error(Limit,Plan,PlanCost,best_plan_unbounded).    
+
+%%%
+is_tabled_state(S) =>
+    bp.b_IS_PLANNER_STATE_c(S).
+
+throw_plan_arg_error(_Limit,Plan,_PlanCost,Source), nonvar(Plan) =>
+    handle_exception($var_expected(Plan), Source).
+throw_plan_arg_error(_Limit,_Plan,PlanCost,Source), nonvar(PlanCost) =>
+    handle_exception($var_expected(PlanCost), Source).
+throw_plan_arg_error(Limit,_Plan,_PlanCost,Source), integer(Limit) =>
+    handle_exception($nonnegative_integer_expected(Limit), Source).
+throw_plan_arg_error(Limit,_Plan,_PlanCost,Source) =>
+    handle_exception($integer_expected(Limit), Source).
+
diff --git a/picat/sys.pi b/picat/sys.pi
new file mode 100644
index 0000000..a336094
--- /dev/null
+++ b/picat/sys.pi
@@ -0,0 +1,35 @@
+module sys.  % imported by default
+
+abort => abort.
+cl(File) => cl(File).
+cl_facts(Facts) => cl_facts(Facts,[]).
+cl_facts(Facts, IndexInfo) => cl_facts(Facts,IndexInfo).
+cl_facts_table(Facts) => cl_facts_table(Facts,[]).
+cl_facts_table(Facts, IndexInfo) => cl_facts_table(Facts,IndexInfo).
+command(String) = command(String).
+compile(File) => compile(File).
+compile_bp(File) => compile_bp(File).
+compile_files_to_c(Fs, CFile) => compile_files_to_c(Fs,CFile).
+debug => debug.
+exit => halt.
+garbage_collect => garbage_collect.
+garbage_collect(Size) => garbage_collect(Size).
+halt => halt.
+help => help.
+initialize_table => initialize_table.
+load(File) => load(File).
+loaded_modules() = loaded_modules().
+nodebug => nodebug.
+nolog => nolog.
+nospy => nospy.
+notrace => notrace.
+picat_path() = picat_path().
+spy(Point) => spy(Point).
+statistics => statistics.
+statistics(Name, Val) => statistics(Name,Val).
+statistics_all() = statistics_all().
+time(Goal) => throw($meta_meta_call_not_allowed(time(Goal))).
+time2(Goal) => throw($meta_meta_call_not_allowed(time2(Goal))).
+time_out(Goal,Limit, Res) => throw($meta_meta_call_not_allowed(time_out(Goal,Limit,Res))).
+trace => trace.
+
diff --git a/picat/util.pi b/picat/util.pi
new file mode 100644
index 0000000..6d6b70e
--- /dev/null
+++ b/picat/util.pi
@@ -0,0 +1,294 @@
+/*
+  Some general utilities in Picat.
+  Others have been moved to the modules basic, math, sys, etc.
+  by Hakan Kjellerstrand and Neng-Fa Zhou.
+*/
+
+module util.
+
+% array_matrix_to_list(A) = array_matrix_to_list(A).
+% array_matrix_to_list_matrix(A) = array_matrix_to_list_matrix(A).
+% chunks_of(L,N) = chunks_of(L,N).
+% columns(Matrix) = columns(Matrix).
+% diagonal1(Matrix) = diagonal1(Matrix).
+% diagonal2(Matrix) = diagonal2(Matrix).
+% drop(L,N) = drop(L,N).
+% find(String, SubString, From, To) => find(String, SubString, From, To).
+% find_first_of(Compound,Pattern) = find_first_of(Compound,Pattern).
+% find_ignore_case(String, SubString, From, To) => find_ignore_case(String, SubString, From, To).
+% find_last_of(Compound,Pattern) = find_last_of(Compound,Pattern).
+% join(S) = join(S).
+% join(S, Seperator) = join(S, Seperator).
+% list_matrix_to_array_matrix(L) = list_matrix_to_array_matrix(L).
+% lstrip(L) = lstrip(L," \t\n\r").
+% lstrip(L,Chars) = lstrip(L,Chars).
+% matrix_multi(A,B) = matrix_multi(A,B).
+% nextto(X,Y,List) => nextto(X,Y,List).
+% permutation(Xs,Ys) => permutation(Xs,Ys).
+% permutations(Xs) = permutations(Xs).
+% power_set(Set) = power_set(Set).
+% replace(Term,Old,New) = replace(Term,Old,New).
+% replace_at(Compound,I,NewArg) = replace_at(Compound,I,NewArg).
+% rows(Matrix) = rows(Matrix).
+% rstrip(L) = rstrip(L," \t\n\r").
+% rstrip(L,Chars) = rstrip(L,Chars).
+% split(Str) = split(Str).
+% split(Str,Seperators) = split(Str,Seperators).
+% strip(L) = strip(L," \t\n\r").
+% strip(L,Chars) = strip(L,Chars).
+% take(L,N) = take(L,N).
+% transpose(Matrix) = transpose(Matrix).
+
+%
+% Convert a 2D array to a list
+%
+array_matrix_to_list(A) = L =>
+    NRows = A.length,
+    NCols = A[1].length,
+    L = [A[I,J] : I in 1..NRows, J in 1..NCols].
+
+% Convert a 2D array to a 2D matrix of lists
+array_matrix_to_list_matrix(A) = L =>
+    L = [A[I].to_list() : I in 1..A.length].
+
+% Convert a 2D list matrix to a 2D array matrix
+list_matrix_to_array_matrix(L) = A =>
+    bp.list_matrix_to_array_matrix(L,A).
+
+%
+% Join a list of strings with a join character.
+% Res = join(String,JoinChar)
+%
+join(S) = join(S," ").
+
+join(S, JoinAtm) = Res, atom(JoinAtm) => 
+    join_aux(S, JoinAtm.to_string(), Res).
+join(S,JoinStr) = Res, string(JoinStr) => 
+    join_aux(S,JoinStr,Res).
+join(S,JoinStr) = _Res => 
+    handle_exception($string_or_atom_expected(JoinStr), $join(S,JoinStr)).
+
+private
+join_aux([],_JoinStr,Res)  => Res = [].
+join_aux([W|Str],JoinStr,Res)  => 
+    once(append(W,Res1,Res)),
+    (Str == [] -> 
+    Res2 = Res1
+    ;
+        once(append(JoinStr,Res2,Res1))
+    ),
+    join_aux(Str,JoinStr,Res2).
+
+%%%%
+% replace occurrences of Old in T (a variable or an atomic value) by New
+replace(T,Old,New) = Res =>
+    replace_aux(T,Old,New,Res).
+    
+replace_aux(Old,Old,New,Res) => Res = New.
+replace_aux(T,Old,New,Res), atomic(T) => Res = T.
+replace_aux(T,Old,New,Res), var(T) => Res = T.    
+replace_aux([H|T],Old,New,Res) =>
+    Res = [NH|NT],
+    replace_aux(H,Old,New,NH), 
+    replace_aux(T,Old,New,NT).
+replace_aux(T,Old,New,Res) =>
+    Res = new_struct(T.name,T.length),
+    foreach(I in 1 .. T.length) 
+        replace_aux(T[I],Old,New,Res[I]) 
+    end.
+
+%%%%
+% return a copy of the compound value, replacing the Ith argument by NewVal
+replace_at(List,I,NewVal) = NewList, integer(I), list(List) =>
+    replace_list_at(List,I,NewVal,NewList,ErrorFlag),
+    (var(ErrorFlag) -> true; handle_exception($domain_error(I), $replace_at)).
+replace_at(Struct,I,NewVal) = NewStruct, integer(I), struct(Struct) =>
+    Arity = len(Struct),
+    (I >= 1, I =< Arity ->
+    NewStruct = new_struct(Struct.name,Arity),
+    foreach(J in 1 .. Arity) 
+            (J == I ->
+            NewStruct[J] = NewVal
+        ;
+            NewStruct[J] = Struct[J]
+        )
+    end
+    ;
+        handle_exception($domain_error(I), $replace_at)
+    ).
+
+private
+replace_list_at([_|List],1,NewVal,NewList,_ErrorFlag) => NewList = [NewVal|List].
+replace_list_at([E|List],I,NewVal,NewList,ErrorFlag) =>
+    NewList = [E|NewList1],
+    replace_list_at(List,I-1,NewVal,NewList1,ErrorFlag).
+replace_list_at(_,_,_,_,ErrorFlag) => ErrorFlag = 1.
+
+% match a string
+%   find(String, SubString,From,To)
+%
+% (If we want to have multiple results it must be a predicate,
+%  not a function.)
+find(String, SubString, From, To) =>
+    SubLen = SubString.length,
+    bp.append(Pre,SubString,_,String),
+    From = Pre.length+1,
+    To = From+SubLen-1.
+
+% Case insensitive match
+find_ignore_case(String, SubString, From, To) =>
+    String2 = String.to_lowercase(),
+    SubString2 = SubString.to_lowercase(),
+    find(String2,SubString2,From,To).
+
+%%%%
+% searches for the first argument that unifies with Pattern and returns the argument's index
+find_first_of(List,Pattern) = Index, list(List) =>
+    find_list_first_of(List,Pattern,1,Index).
+find_first_of(Struct,Pattern) = Index, struct(Struct) =>
+    find_struct_first_of(Struct,Pattern,1, len(Struct), Index).
+find_first_of(Struct,_Pattern) = _ =>
+    handle_exception($compound_expected(Struct), find_first_of).
+
+find_list_first_of([],_Pattern,_CurIndex,Index) => Index = -1.
+find_list_first_of([E|L],Pattern,CurIndex,Index), E != Pattern =>  % not unifiable
+    find_list_first_of(L,Pattern,CurIndex+1,Index).
+find_list_first_of(_L,_Pattern,CurIndex,Index) => Index = CurIndex.
+
+find_struct_first_of(_Struct,_Pattern,CurIndex,Arity,Index), CurIndex > Arity => Index = -1.
+find_struct_first_of(Struct,Pattern,CurIndex,Arity,Index), Struct[CurIndex] != Pattern =>  % not unifiable
+    find_struct_first_of(Struct,Pattern,CurIndex+1,Arity,Index).
+find_struct_first_of(_,_,CurIndex,_Arity,Index) => Index = CurIndex.
+
+%%%%
+% searches for the last argument that unifies with Pattern and returns the argument's index
+find_last_of(List,Pattern) = Index, list(List) =>
+    find_list_last_of(List,Pattern,1,-1,Index).
+find_last_of(Struct,Pattern) = Index, struct(Struct) =>
+    find_struct_last_of(Struct,Pattern, len(Struct), Index).
+find_last_of(Struct,_) = _ =>
+    handle_exception($compound_expected(Struct), find_last_of).
+
+find_list_last_of([],_Pattern,_CurIndex,Index0,Index) => Index = Index0.
+find_list_last_of([E|L],Pattern,CurIndex,Index0,Index), E != Pattern =>  % not unifiable
+    find_list_last_of(L,Pattern,CurIndex+1,Index0,Index).
+find_list_last_of([_|L],Pattern,CurIndex,_Index0,Index) =>
+    find_list_last_of(L,Pattern,CurIndex+1,CurIndex,Index).
+
+find_struct_last_of(_Struct,_Pattern,0,Index) => Index = -1.
+find_struct_last_of(Struct,Pattern,CurIndex,Index), Struct[CurIndex] != Pattern =>  % not unifiable
+    find_struct_last_of(Struct,Pattern,CurIndex-1,Index).
+find_struct_last_of(_,_,CurIndex,Index) => Index = CurIndex.
+
+% A*B=C
+matrix_multi(A,B) = C, array(A), array(B) =>   % A and B must be array matricies
+    C = new_array(A.length,B[1].length),
+    foreach(I in 1..A.length, J in 1..B[1].length)
+        C[I, J] = sum([A[I,K]*B[K,J] : K in 1..A[1].length])
+    end.
+
+% nondet
+permutation(Xs, Ys) => bp.permutation(Xs,Ys).
+
+% generate permutations
+permutations([]) = [[]].
+permutations([H|T]) = [insert(P,I,H) : P in Ps, I in 1..P.length+1] => Ps = permutations(T).
+
+%nonet
+nextto(X,Y, List) => bp.nextto(X,Y,List).
+
+% generate the power set
+power_set([]) = [[]].
+power_set([H|T]) = P1++P2 =>
+    P1 = power_set(T),
+    P2 = [[H|S] : S in P1].
+
+%
+% Split a string into tokens given some split chars
+% List = split(String, Seperators)
+%
+split(Str) = split(Str," \t\n\r").   % use white spaces as the default set of separators
+
+split(Str,Seperators) = Tokens =>
+    bp.picat_split_string(Str,Seperators,Tokens).
+
+lstrip(L) = lstrip(L," \t\n\r").
+
+lstrip([],_Elms) = [].
+lstrip([E|L],Elms) = NewL, membchk(E, Elms) => NewL = lstrip(L,Elms).
+lstrip(L,_Elms) = NewL => NewL = L.
+
+rstrip(L) = rstrip(L," \t\n\r").
+
+rstrip(L, Elms) = L.reverse().lstrip(Elms).reverse().
+
+strip(L) = strip(L," \t\n\r").
+
+strip(L, Elms) = L.lstrip(Elms).rstrip(Elms).
+
+%
+% Transpose a 2D matrix
+%
+transpose(Matrix) = Transposed, array(Matrix) =>   % array matrix
+    N = Matrix.length,
+    M = Matrix[1].length,
+    Transposed = new_array(M,N),
+    foreach(I in 1..N, J in 1..M)
+        Transposed[J,I] = Matrix[I,J]
+    end.
+transpose(Matrix) = Transposed =>   % assumed to be list matrix
+    N = Matrix.length,
+    M = Matrix[1].length,
+    Transposed = [Mj : J in 1..M, Mj = [Matrix[I,J] : I in 1..N]].
+
+%=============
+% for matrices (inspired by B-Prolog's ^rows, ^columns, ^diag1, ^diag2)
+%
+% These should be put in util.pi since transpose/1 is used.
+
+rows(M) = M, list(M)  => true.
+rows(A) = Rows, array(A)  =>  NRows = A.length, Rows = [A[I] : I in 1..NRows].
+
+columns(M) = [Column : Column in M.transpose()], list(M) => true.
+columns(A) = [Column : Column in A.transpose()], array(A) => true.
+
+diagonal1(M) = [M[I,I] : I in 1..M.length], list(M)  => true.
+diagonal1(A) = [A[I,I] : I in 1..A.length], array(A)  => true.
+
+diagonal2(M) = [M[I,M.length-I+1] : I in 1..M.length], list(M)  => true.
+diagonal2(A) = [A[I,A.length-I+1] : I in 1..A.length], array(A)  => true.
+
+
+%=============
+% from Haskell prelude
+%
+take(L,N) = Taken, list(L), integer(N), take_aux(L,N,Taken) => true.
+
+private
+take_aux([H|T],N,Taken), N > 0 => Taken = [H|TakenR], take_aux(T,N-1,TakenR).
+take_aux(_List,_N,Taken) => Taken = [].
+
+drop(L,N) = Taken, list(L), integer(N), drop_aux(L,N,Taken) => true.
+
+private
+drop_aux([_|T],N,Taken), N > 0 => drop_aux(T,N-1,Taken).
+drop_aux(L,_,Taken), list(L) => Taken = L.
+
+chunks_of([],_N) = [].
+chunks_of(L,N) = Chunks, list(L) =>
+    Chunks = [Chunk|ChunksR],
+    chunks_of(L,Chunk,0,N,ChunksR).
+
+chunks_of([],Chunk,_,_,Chunks) => Chunk = [], Chunks = [].
+chunks_of(L,Chunk,N,N,Chunks) =>
+    Chunk = [],
+    (L == [] ->
+        Chunks = []
+    ;
+        Chunks = [NextChunk|ChunksR],
+        chunks_of(L,NextChunk,0,N,ChunksR)
+    ).
+chunks_of([X|Xs],Chunk,Count,N,Chunks) =>
+    Chunk = [X|ChunkR],
+    chunks_of(Xs,ChunkR,Count+1,N,Chunks).
+