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

`is_element/Q` := (A::set) -> proc(t)
 local A3,abc,a,b,c,r2;
 global reason;

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

 A3 := `list_elements/triples`(A);

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

 for abc in A3 do
  a,b,c := op(abc);
  if not(`is_element/RR`(t[a,b,c])) then
   reason := [convert(procname,string),"t[a,b,c] is not in R",a,b,c,t[a,b,c]];
   return false;
  fi;
  if not(is(t[a,b,c] >= 0)) then return false; fi;
 od;

 for abc in A3 do
  a,b,c := op(abc);
  if simplify(t[a,b,c] - t[b,a,c]) <> 0 then
   reason := [convert(procname,string),"t[a,b,c] <> t[b,a,c]",a,b,c,t[a,b,c],t[b,a,c]];
   return false;
  fi;
  r2 := simplify(t[a,b,c]^2 + t[b,c,a]^2 + t[c,a,b]^2);
  if simplify(r2-1) <> 0 then
   reason := [convert(procname,string),"t is not normalised on {a,b,c}",a,b,c,r2];
   return false;
  fi;
 od;

 return true;
end:

`psi/Fbar/Q` := (N::posint) -> (A::set) -> proc(x)
 local A3,t,abc,a,b,c,f,g;

 A3 := `list_elements/triples`(A);
 t := table();

 for abc in A3 do
  a,b,c := op(abc);
  f := `normalise_2/F`(N)({a,b,c})(x[{a,b,c}]);
  g := `normalise_2/F`(N)({a,b})(x[{a,b}]);
  t[a,b,c] := combine(simplify(`dot/R`(N)(f[b] -~ f[a],g[b])*sqrt(2/3)));
 od;

 return(eval(t));
end:

ext: mpl


text:
######################################################################
# Fbar(N)(A) is the Fulton-MacPherson compactification of F(N)(A),
# defined in terms of coherent families as in the LaTeX document.

`is_element/Fbar` := (N::posint) -> (A::set) -> proc(x)
 local B,C,P,Q;
 global reason;

 P := {op(`list_elements/big_subsets`(A))};

 if not(is_table_on(P)(x)) then
  reason := [convert(procname,string),"x is not a table indexed by big subsets of A",eval(x)];
  return false;
 fi;

 for B in P do
  if not(`is_element/SW`(N)(B)(x[B])) then
   reason := [convert(procname,string),"x[B] is not in SW(N)(B)",x[B],N,B,reason];
   return false;
  fi;

  Q := {op(`list_elements/big_subsets`(B))} minus {B};

  for C in Q do
   if not(`is_element/SWW`(N)(B,C)([x[B],x[C]])) then
    reason := [convert(procname,string),"(x[B],x[C]) is not in SWW(N)(B,C)",x[B],x[C],N,B,C,reason];
    return false;
   fi;
  od;
 od;

 return true;
end;

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

`normalise/Fbar` := (N::posint) -> (A::set) -> proc(x)
 local x0,B,P;

 x0 := table();
 P := `list_elements/big_subsets`(A);

 for B in P do x0[B] := `normalise/SW`(N)(B)(x[B]); od;
 return eval(x0);
end;

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

`bottom_normalise/Fbar` := (N::posint) -> (A::set) -> proc(x)
 local x0,B,P;

 x0 := table();
 P := `list_elements/big_subsets`(A);

 for B in P do x0[B] := `bottom_normalise/SW`(N)(B)(x[B]); od;
 return eval(x0);
end;

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

`is_equal/Fbar` := (N::posint) -> (A::set) -> proc(x,y)
 local B,P;
 global reason;

 P := `list_elements/big_subsets`(A);

 for B in P do
  if not(`is_equal/SW`(N)(B)(x[B],y[B])) then
   reason := [convert(procname,string),"x[B] <> y[B]",B,x[B],y[B],reason];
   return false;
  fi;
 od;

 return true;
end;

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

`is_leq/Fbar` := NULL;

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

`is_interior/Fbar` := (N::posint) -> (A::set) -> proc(x)
 local V;

 V := map(a -> x[A][a],A);
 return evalb(nops(V) = nops(A));
end;

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

`inc/F/Fbar` := (N::posint) -> (A::set) -> proc(u)
 local x,P,T,t;

 x := table();
 P := `list_elements/big_subsets`(A);
 for T in P do
  x[T] := table();
  for t in T do
   x[T][t] := u[t];
  od;
 od;

 return eval(x);
end:

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

`res/Fbar/F` := (N::posint) -> (A::set) -> proc(x)
 if not(`is_interior/Fbar`(N)(A)(x)) then
  return FAIL;
 fi;

 return x[A];
end;

######################################################################
# The functions below could be made much more efficient, but for
# the moment we will stick with the most obvious approach.

`is_critical/Fbar` := (N::posint) -> (A::set) -> proc(x,U)
 local T,TT;

 if nops(U) = 1 then return true; fi;

 TT := map(E -> sort(U union E),
           `list_elements/nonempty_subsets`(A minus U));

 for T in TT do
  if not(`is_zero/W`(N)(U)(`rho/W`(T,U)(x[T]))) then
   return false;
  fi;
 od;

 return true;
end;

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

`critical_tree/Fbar` := (N::posint) -> (A::set) -> proc(x)
 select(U -> `is_critical/Fbar`(N)(A)(x,U),
        {op(`list_elements/nonempty_subsets`(A))});
end;

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

`random_element/Fbar` := (N::posint) -> (A::set) -> proc()
 local TT,TS,parent,children,u,v,w,t,y,P,T,U,C;

 TT := `random_element/trees`(A)();
 TS := sort([op(TT)],(A,B) -> (nops(A) < nops(B)));
 parent := `parent_map`(A)(TT);
 children := `children_map`(A)(TT);

 u := table();
 for T in TS do
  u[T] := `random_element/F`(N)(children[T])();
  v[T] := table();
  for C in children[T] do
   y := u[T][C];
   for t in C do
    v[T][t] := y;
   od;
  od;
 od:

 w := table();

 P := {op(`list_elements/big_subsets`(A))};

 for U in P do
  for T in TT do
   if U minus T = {} then
    w[U] := table();
    for t in U do
     w[U][t] := v[T][t];
    od;
    break;
   fi;
  od;
 od;

 return eval(w);
end;

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

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

`eta/Fbar` := (N::posint) -> (A::set) -> `if`(nops(A) = 1,table(),FAIL);

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

`gamma/Fbar` := (N::posint) -> (A::set,B::set) -> (p) -> proc(y,x)
 local z,TT,T,U,m,a;

 z := table();
 TT := `list_elements/big_subsets`(A);
 for T in TT do
  U := map(a -> p[a],T);
  if nops(U) > 1 then
   m := y[U];
   z[T] := table({seq(a = m[p[a]],a in T)});
  else
   z[T] := eval(x[op(U)][T]);
  fi;
 od;

 return eval(z);
end;


text:
TITLE
    ShortPath;

INDEX
    node := DATAFILE("Cities.dat", 1);

    FromNode := node;
    ToNode   := node;

DATA
    Distance[FromNode, ToNode] :=  SPARSEFILE("Routes.dat");

    StartCity  :=  "Washington";
    EndCity    :=  "San Francisco";
    Flow       :=  1.0;

VARIABLES
    Path[FromNode, ToNode] WHERE (Distance > 0);

MODEL

    MIN TotalDist = SUM(FromNode, ToNode: Distance * Path);

SUBJECT TO

    FlowBalance[node] :

      Flow IF (node=StartCity) - Flow IF (node=EndCity) +

       SUM(FromNode: Path[FromNode, ToNode:=node])
      =
       SUM(ToNode: Path[FromNode:=node, ToNode]);

END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


text:
(*
 Copyright (C) 2021 Susi Lehtola

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

(* type: mgga_exc *)
$include "mgga_c_scan.mpl"

(* interpolation function, defined in main text *)
rregtm_gamma1 := 0.2:
rregtm_g := a -> (1 + rregtm_gamma1)*a / (rregtm_gamma1 + a):
rregtm_f2g := g -> 3*g^3 / (1 + g^3 + g^6):
rregtm_f2 := a -> rregtm_f2g(rregtm_g(a)):

(* equation 3 *)
rregtm_f := (rs, z, xt, xs0, xs1, ts0, ts1) ->
  f_pbe(rs, z, xt, xs0, xs1) + rregtm_f2(scan_alpha(z, xt, ts0, ts1))*(
    + scan_e0(rs, z, X2S*2^(1/3)*xt)
    - f_pbe(rs, z, xt, xs0, xs1)
  ):

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


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

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

(* type: work_gga_c *)

$define lda_c_pz_params
$include "lda_c_pz.mpl"

malpha := 0.023266:
mbeta  := 7.389e-6:
mgamma := 8.723:
mdelta := 0.472:
aa     := 0.001667:
bb     := 0.002568:
ftilde := 1.745*0.11:

(* Equation (4) *)
DD := z  -> sqrt((1 + z)^(5/3) + (1 - z)^(5/3))/sqrt(2):

(* Equation (6) *)
CC := rs ->
  + aa
  + (bb + malpha*rs + mbeta*rs^2)/(1 + mgamma*rs + mdelta*rs^2 + 1.0e4*mbeta*rs^3):
CCinf := aa + bb:

(* Equation (9) *)
x1   := (rs, xt) -> xt/sqrt(rs/RS_FACTOR):
mPhi := (rs, xt) -> ftilde*(CCinf/CC(rs))*x1(rs, xt):

(* Equation (8) *)
H := (rs, z, xt) -> x1(rs, xt)^2*exp(-mPhi(rs, xt))*CC(rs)/DD(z):

f := (rs, z, xt, xs0, xs1) ->
  f_pz(rs, z) + H(rs, z, xt):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms