text:
#############################################################################
##
#W  treeautgrp.gi              automgrp package                Yevgen Muntyan
#W                                                             Dmytro Savchuk
##
#Y  Copyright (C) 2003 - 2018 Yevgen Muntyan, Dmytro Savchuk
##

###############################################################################
##
#M  TreeAutomorphismGroup(<G>, <S>)
##
##  Wreath product
##
InstallOtherMethod(TreeAutomorphismGroup, [IsTreeAutomorphismGroup, IsSymmetricGroup],
function(G, S)
  local g_gens, s_gens, g_deg, s_deg, orbs, o, gens;
  s_deg := Maximum(MovedPointsPerms(S));
  g_deg := DegreeOfTree(G);
  gens := List(GeneratorsOfGroup(S), s -> [List([1..s_deg], i -> One(G)), s]);
  Append(gens, List(GeneratorsOfGroup(G), g ->
                [Concatenation([g], List([2..s_deg], i -> One(g))), ()]));
  Apply(gens, g -> TreeAutomorphism(g[1], g[2]));
  return GroupWithGenerators(gens);
end);

###############################################################################
##
#M  UseSubsetRelation(<super>, <sub>)
##
InstallMethod(UseSubsetRelation, "for [IsTreeAutomorphismGroup, IsTreeAutomorphismGroup]",
              [IsTreeAutomorphismGroup, IsTreeAutomorphismGroup],
function(super, sub)
  if HasIsSphericallyTransitive(super) then
    if not IsSphericallyTransitive(super) then
      Info(InfoAutomGrp, 3, "IsSphericallyTransitive(sub): false");
      Info(InfoAutomGrp, 3, "  super is not spherically transitive");
      Info(InfoAutomGrp, 3, "  super = ", super);
      Info(InfoAutomGrp, 3, "  sub = ", sub);
      SetIsSphericallyTransitive(sub, false); fi; fi;

  TryNextMethod();
end);

###############################################################################
##
#M  SphericalIndex(<G>)
##
InstallMethod(SphericalIndex, "for [IsTreeAutomorphismGroup]",
              [IsTreeAutomorphismGroup],
function(G)
  return SphericalIndex(GeneratorsOfGroup(G)[1]);
end);
InstallMethod(TopDegreeOfTree, "for [IsTreeAutomorphismGroup]",
              [IsTreeAutomorphismGroup],
function(G)
  return TopDegreeOfTree(GeneratorsOfGroup(G)[1]);
end);
InstallMethod(DegreeOfTree, "for [IsTreeAutomorphismGroup]",
              [IsTreeAutomorphismGroup],
function(G)
  return DegreeOfTree(GeneratorsOfGroup(G)[1]);
end);

###############################################################################
##
#M  IsSphericallyTransitive (<G>)
#M  CanEasilyTestSphericalTransitivity (<G>)
##
##  Fractalness implies spherical transitivity.
##
InstallTrueMethod(IsSphericallyTransitive, IsFractal);
InstallTrueMethod(CanEasilyTestSphericalTransitivity, HasIsSphericallyTransitive);

InstallImmediateMethod(IsSphericallyTransitive, IsTreeAutomorphismGroup and HasIsFinite, 0,
function(G)
  if IsFinite(G) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
    Info(InfoAutomGrp, 3, "  G is finite");
    Info(InfoAutomGrp, 3, "  G = ", G);
    return false;
  fi;
  TryNextMethod();
end);

InstallMethod(IsSphericallyTransitive, "for [IsTreeAutomorphismGroup]",
              [IsTreeAutomorphismGroup],
function (G)
  local i, k, stab;

  if DegreeOfTree(G)=1 then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): true");
    Info(InfoAutomGrp, 3, "  G acts on 1-ary tree");
    return true;
  fi;

  if HasIsFinite(G) and IsFinite(G) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
    Info(InfoAutomGrp, 3, "  G is finite");
    Info(InfoAutomGrp, 3, "  G = ", G);
    return false;
  fi;

  if not IsTransitiveOnLevel(G, 1) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
    Info(InfoAutomGrp, 3, "  G is not transitive on ", i, "-th level");
    Info(InfoAutomGrp, 3, "  G = ", G);
    return false;
  fi;

  if IsActingOnBinaryTree(G) then
    if AG_AbelImageSpherTrans() in AbelImage(G) then
      Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): true");
      Info(InfoAutomGrp, 3, "  using AbelImage");
      Info(InfoAutomGrp, 3, "  G = ", G);
      return true;
    fi;
  fi;

  if not IsAutomGroup(G) then
    stab := StabilizerOfVertex(G, 1);
    return IsSphericallyTransitive(Projection(stab, 1));
  fi;

#  TryNextMethod();
  return fail;
end);

###############################################################################
##
#M  IsTransitiveOnLevel (<G>, <k>)
##
InstallMethod(IsTransitiveOnLevel, [IsTreeAutomorphismGroup, IsPosInt],
function(G, k)
  return IsTransitive(G, AG_TreeLevelTuples(SphericalIndex(G), k));
end);

###############################################################################
##
#M  AbelImage (<G>)
##
InstallMethod(AbelImage, [IsTreeAutomorphismGroup],
function(G)
  local ab, gens, abgens;

  abgens := List(GeneratorsOfGroup(G), g -> AbelImage(g));
  return AdditiveGroup(abgens);
end);

###############################################################################
##
#M  IsFractal (<G>)
##
##  Fractalness implies spherical transitivity, hence not spherical transitive
##  group isn't fractal.
##
InstallImmediateMethod(IsFractal, HasIsSphericallyTransitive, 0,
function(G)
  if not IsSphericallyTransitive(G) then return false; fi;
  TryNextMethod();
end);

InstallImmediateMethod(IsFractal, IsTreeAutomorphismGroup and HasIsFinite, 0,
function(G)
  if IsFinite(G) and DegreeOfTree(G)>1 then return false; fi;
  TryNextMethod();
end);

InstallMethod(IsFractal, "for [IsTreeAutomorphismGroup]", [IsTreeAutomorphismGroup],
function(G)
  local i;

  if DegreeOfTree(G)=1 then
    SetIsSphericallyTransitive(G, true);
    Info(InfoAutomGrp, 3, "IsFractal(G): true");
    Info(InfoAutomGrp, 3, "  G acts on 1-ary tree");
    return true;
  fi;

  if not IsTransitive(PermGroupOnLevel(G, 1), [1..DegreeOfTree(G)]) then
    SetIsSphericallyTransitive(G, false);
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
    Info(InfoAutomGrp, 3, "IsFractal(G): false");
    Info(InfoAutomGrp, 3, "  G is not transitive on first level");
    return false;
  fi;

  if HasIsFinite(G) and IsFinite(G) then
    Info(InfoAutomGrp, 3, "IsFractal(G): false");
    Info(InfoAutomGrp, 3, "  G is finite");
    return false;
  fi;

  if CanEasilyTestSphericalTransitivity(G) and not IsSphericallyTransitive(G) then
    Info(InfoAutomGrp, 3, "IsFractal(G): false");
    Info(InfoAutomGrp, 3, "  G is not spherically transitive");
    return false;
  fi;

  if ProjStab(G, 1) <> G then
    Info(InfoAutomGrp, 3, "IsFractal(G): false");
    Info(InfoAutomGrp, 3, "  ProjStab(G, 1) <> G");
    return false;
  fi;
  Info(InfoAutomGrp, 3, "IsFractal(G): true");
  Info(InfoAutomGrp, 3, "  ProjStab(G, 1) = G and G is transitive on first level");
  return true;
end);

###############################################################################
##
#M  Size (<G>)
##
InstallImmediateMethod(Size, HasIsFractal, 0,
function(G)
  if IsFractal(G) and DegreeOfTree(G)>1 then return infinity; fi;
  TryNextMethod();
end);

InstallImmediateMethod(Size, HasIsSphericallyTransitive, 0,
function(G)
  if IsSphericallyTransitive(G) and DegreeOfTree(G)>1 then return infinity; fi;
  TryNextMethod();
end);

###############################################################################
##
#M  PermGroupOnLevelOp (<G>, <k>)
##
InstallMethod(PermGroupOnLevelOp, "for [IsTreeAutomorphismGroup, IsPosInt]",
              [IsTreeAutomorphismGroup, IsPosInt],
function(G, k)
  local pgens, pgroup;

  pgens := List(GeneratorsOfGroup(G), g -> PermOnLevel(g, k));
  if pgens = [] then
    return Group(());
  else
    pgroup := GroupWithGenerators(pgens);
    if IsActingOnBinaryTree(G) then
      SetIsPGroup(pgroup, true);
      if IsTrivial(pgroup) then
        SetPrimePGroup(pgroup, fail);
      else
        SetPrimePGroup(pgroup, 2);
      fi;
    fi;
    return pgroup;
  fi;
end);

###############################################################################
##
#M  StabilizerOfFirstLevel(G)
##
InstallMethod(StabilizerOfFirstLevel,
              "for [IsTreeAutomorphismGroup]",
              [IsTreeAutomorphismGroup],
function (G)
  return StabilizerOfLevel(G, 1);
end);

###############################################################################
##
#M  StabilizerOfLevelOp(G, k)
##
InstallMethod(StabilizerOfLevelOp,
              "for [IsTreeAutomorphismGroup, IsPosInt]",
              [IsTreeAutomorphismGroup, IsPosInt],
function (G, k)
  local perms, S, F, hom, chooser, s, f, gens;

  if DegreeOfTree(G)=1 then
    return G;
  fi;

  if FixesLevel(G, k) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
    Info(InfoAutomGrp, 3, "  G is not transitive on level", k);
    Info(InfoAutomGrp, 3, "  G = ", G);
    SetIsSphericallyTransitive(G, false);
    return G;
  fi;

  perms := List(GeneratorsOfGroup(G), a -> PermOnLevel(a, k));
  S := GroupWithGenerators(perms);
  F := FreeGroup(Length(perms));
  hom := GroupHomomorphismByImagesNC(F, S, GeneratorsOfGroup(F), perms);
  gens := FreeGeneratorsOfGroup(Kernel(hom));
  gens := List(gens, w  ->
    MappedWord(w, GeneratorsOfGroup(F), GeneratorsOfGroup(G)));
  gens := __AG_SimplifyGroupGenerators(gens);
  if IsEmpty(gens) then
    return TrivialSubgroup(G);
  else
    return SubgroupNC(G, gens);
  fi;
end);

###############################################################################
##
#M  StabilizerOfVertex(G, k)
##
InstallMethod(StabilizerOfVertex,
              "for [IsTreeAutomorphismGroup, IsPosInt]",
              [IsTreeAutomorphismGroup, IsPosInt],
function (G, k)
  local X, S, F, hom, s, f, gens, i, action;

  if k > TopDegreeOfTree(G) then
    Print("error in StabilizerOfVertex(IsTreeAutomorphismGroup, IsPosInt):\n");
    Print("  given k is not a valid vertex\n");
    return fail;
  fi;
  if FixesVertex(G, k) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
    Info(InfoAutomGrp, 3, "  G fixes vertex", k);
    Info(InfoAutomGrp, 3, "  G = ", G);
    Info(InfoAutomGrp, 3, "  k = ", k);
    SetIsSphericallyTransitive(G, false);
    return G;
  fi;
  if TopDegreeOfTree(G) = 2 then
    return StabilizerOfFirstLevel(G);
  fi;

  X := List(GeneratorsOfGroup(G), a -> Perm(a));
  S := GroupWithGenerators(X);
  F := FreeGroup(Length(X));
  hom := GroupHomomorphismByImagesNC(F, S,
              GeneratorsOfGroup(F), X);
  action := function(k, w)
      return k^Image(hom, w);
  end;
  gens := FreeGeneratorsOfGroup(Stabilizer(F, k, action));
  gens := List(gens, w  ->
    MappedWord(w, GeneratorsOfGroup(F), GeneratorsOfGroup(G)));
  gens := __AG_SimplifyGroupGenerators(gens);
  if IsEmpty(gens) then
    return TrivialSubgroup(G);
  else
    return SubgroupNC(G, gens);
  fi;
end);

###############################################################################
##
#M  StabilizerOfVertex(G, seq)
##
InstallMethod(StabilizerOfVertex,
              "for [IsTreeAutomorphismGroup, IsList]",
              [IsTreeAutomorphismGroup, IsList],
function (G, seq)
  local X, S, F, hom, s, f, gens, stab, action, i, v;

  if Length(seq) = 0 then
    return G;
  fi;

  if TopDegreeOfTree(G) = 2 and Length(seq) = 1 then
    return StabilizerOfFirstLevel(G);
  fi;

  if FixesVertex(G, seq) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
    Info(InfoAutomGrp, 3, "  G fixes vertex", seq);
    Info(InfoAutomGrp, 3, "  G = ", G);
    SetIsSphericallyTransitive(G, false);
    return G;
  fi;

  v := Position(AsList(Tuples([1..DegreeOfTree(G)], Length(seq))), seq);

  X := List(GeneratorsOfGroup(G), a -> [Word(a), PermOnLevel(a, Length(seq))]);
  S := Group(List(X, x -> x[2]));
  F := FreeGroup(Length(X));
  hom := GroupHomomorphismByImagesNC(F, S,
              GeneratorsOfGroup(F), List(X, x -> x[2]));
  action := function(k, w)
      return k^Image(hom, w);
  end;
  gens := FreeGeneratorsOfGroup(Stabilizer(F, v, action));
  gens := List(gens, w  ->
    MappedWord(w, GeneratorsOfGroup(F), GeneratorsOfGroup(G)));
  gens := __AG_SimplifyGroupGenerators(gens);
  if IsEmpty(gens) then
    return TrivialSubgroup(G);
  else
    return SubgroupNC(G, gens);
  fi;
end);

###############################################################################
##
#M  FixesLevel(<G>, <k>)
##
InstallMethod(FixesLevel, "for [IsTreeAutomorphismGroup, IsPosInt]",
              [IsTreeAutomorphismGroup, IsPosInt],
function(G, k)
  if IsTrivial(PermGroupOnLevel(G, k)) then
    if DegreeOfTree(G)>1 then
      Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
      SetIsSphericallyTransitive(G, false);
    fi;
    Info(InfoAutomGrp, 3, "  G fixes level", k);
    Info(InfoAutomGrp, 3, "  G = ", G);
    return true;
  else
    return false;
  fi;
end);

###############################################################################
##
#M  FixesVertex(<G>, <v>)
##
InstallMethod(FixesVertex, "for [IsTreeAutomorphismGroup, IsPosInt]",
              [IsTreeAutomorphismGroup, IsObject],
function(G, v)
  local gens, g;
  gens := GeneratorsOfGroup(G);
  for g in gens do
    if not FixesVertex(g, v) then return false; fi;
  od;
  if DegreeOfTree(G)>1 then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(G): false");
    Info(InfoAutomGrp, 3, "  G fixes vertex", v);
    Info(InfoAutomGrp, 3, "  G = ", G);
    SetIsSphericallyTransitive(G, false);
  fi;
  return true;
end);

###############################################################################
##
#M  ProjectionOp (<G>, <k>)
##
InstallMethod(ProjectionOp, "for [IsTreeAutomorphismGroup, IsPosInt]",
              [IsTreeAutomorphismGroup, IsPosInt],
function(G, k)
  return Projection(G, [k]);
end);

# TODO: check whether gens are from the same overgroup;
# check degree of tree and stuff
InstallMethod(__AG_SubgroupOnLevel, [IsTreeAutomorphismGroup,
                                    IsList and IsTreeAutomorphismCollection,
                                    IsPosInt],
function(G, gens, level)
  local a;
  if IsEmpty(gens) then
    a := Section(One(G), List([1..level], i -> 1));
    gens := [a];
  fi;
  return Group(gens);
end);

InstallOtherMethod(__AG_SubgroupOnLevel, [IsTreeAutomorphismGroup,
                                    IsList and IsEmpty,
                                    IsPosInt],
function(G, gens, level)
  return Group(Section(One(G), List([1..level], i -> 1)));
end);

InstallMethod(__AG_SimplifyGroupGenerators, "for [IsList and IsTreeAutomorphismCollection]",
                              [IsList and IsTreeAutomorphismCollection],
function(gens)
  if IsEmpty(gens) then
    return [];
  else
    return Difference(gens, [One(gens[1])]);
  fi;
end);

###############################################################################
##
#M  ProjectionNC(<G>, <v>)
##
InstallMethod(ProjectionNC, "for [IsTreeAutomorphismGroup, IsList]",
              [IsTreeAutomorphismGroup, IsList],
function(G, v)
  local gens, pgens, a;

  if IsEmpty(v) then
    return G;
  fi;

  gens := GeneratorsOfGroup(G);
  pgens := List(gens, g -> Section(g, v));
  pgens := __AG_SimplifyGroupGenerators(pgens);

  return __AG_SubgroupOnLevel(G, pgens, Length(v));
end);

###############################################################################
##
#M  Projection (<G>, <v>)
##
InstallOtherMethod(Projection, "for [IsTreeAutomorphismGroup, IsList]",
                   [IsTreeAutomorphismGroup, IsList],
function(G, v)
  if not FixesVertex(G, v) then
    Error("in Projection(G, v): G does not fix vertex v");
  fi;
  return ProjectionNC(G, v);
end);

###############################################################################
##
#M  ProjStab(<G>, <v>)
##
InstallMethod(ProjStab, "for [IsTreeAutomorphismGroup, IsObject]",
              [IsTreeAutomorphismGroup, IsObject],
function(G, v)
  if HasIsFractal(G) and IsFractal(G) then
    return G;
  else
    if IsPosInt(v) then
      v := [v];
    fi;
    return ProjectionNC(StabilizerOfVertex(G, v), v);
  fi;
end);

###############################################################################
##
#M  \= (<G>, <H>)
##
InstallMethod(\=, "for [IsTreeAutomorphismGroup, IsTreeAutomorphismGroup]",
              [IsTreeAutomorphismGroup, IsTreeAutomorphismGroup],
function(G, H)
  if HasIsFinite(G) and HasIsFinite(H) and
    IsFinite(G) <> IsFinite(H) then
      Info(InfoAutomGrp, 3, "\=(IsTreeAutomorphismGroup, IsTreeAutomorphismGroup): false");
      Info(InfoAutomGrp, 3, "  IsFinite(G)<>IsFinite(H)");
      return false;
  fi;

  if HasSize(G) and HasSize(H) and
    Size(G) <> Size(H) then
      Info(InfoAutomGrp, 3, "\=(IsTreeAutomorphismGroup, IsTreeAutomorphismGroup): false");
      Info(InfoAutomGrp, 3, "  Size(G)<>Size(H)");
      return false;
  fi;

  if CanEasilyTestSphericalTransitivity(G) and
    CanEasilyTestSphericalTransitivity(H) and
      IsSphericallyTransitive(G) <> IsSphericallyTransitive(H) then
        Info(InfoAutomGrp, 3, "\=(IsTreeAutomorphismGroup, IsTreeAutomorphismGroup): false");
        Info(InfoAutomGrp, 3, "  IsSphericallyTransitive(G) <> IsSphericallyTransitive(H)");
        return false;
  fi;

  if HasIsFractal(G) and HasIsFractal(H) and
     IsFractal(G) <> IsFractal(H)
  then
    Info(InfoAutomGrp, 3, "\=(IsTreeAutomorphismGroup, IsTreeAutomorphismGroup): false");
    Info(InfoAutomGrp, 3, "  IsFractal(G) <> IsFractal(H)");
    return false;
  fi;

  return IsSubgroup(G, H) and IsSubgroup(H, G);
end);

###############################################################################
##
#M  IsSubset (<G>, <H>)
##
InstallMethod(IsSubset, "for [IsTreeAutomorphismGroup, IsTreeAutomorphismGroup]",
              [IsTreeAutomorphismGroup, IsTreeAutomorphismGroup],
function(G, H)
  local h, gens;

  if IsTrivial(H) then
    Info(InfoAutomGrp, 3, "IsSubset(G, H): true");
    Info(InfoAutomGrp, 3, "  IsTrivial(H)");
    return true;
  elif IsTrivial(G) then
    Info(InfoAutomGrp, 3, "IsSubset(G, H): false");
    Info(InfoAutomGrp, 3, "  not IsTrivial(H) and IsTrivial(G)");
    return false;
  fi;

  if HasIsFinite(G) and HasIsFinite(H) and
    IsFinite(G) and not IsFinite(H) then
      Info(InfoAutomGrp, 3, "IsSubset(G, H): false");
      Info(InfoAutomGrp, 3, "  IsFinite(G) and not IsFinite(H)");
      return false;
  fi;

  if HasSize(G) and HasSize(H) and
    Size(G) < Size(H) then
      Info(InfoAutomGrp, 3, "IsSubset(G, H): false");
      Info(InfoAutomGrp, 3, "  Size(G) < Size(H)");
      return false;
  fi;

  if  CanEasilyTestSphericalTransitivity(G) and
      CanEasilyTestSphericalTransitivity(H) and
      not IsSphericallyTransitive(G) and
      IsSphericallyTransitive(H)
  then
    Info(InfoAutomGrp, 3, "IsSubset(G, H): false");
    Info(InfoAutomGrp, 3, "  not IsSphericallyTransitive(G) and IsSphericallyTransitive(H)");
    return false;
  fi;

  G := AG_ReducedForm(G);
  H := AG_ReducedForm(H);

  gens := GeneratorsOfGroup(H);
  for h in gens do
    if not h in G then
      Info(InfoAutomGrp, 3, "IsSubset(G, H): false");
      Info(InfoAutomGrp, 3, "  not h in G");
      return false;
    fi;
  od;
  return true;
end);

###############################################################################
##
#M  <g> in <G>
##
InstallMethod(\in, "for [IsTreeAutomorphism, IsTreeAutomorphismGroup]",
              [IsTreeAutomorphism, IsTreeAutomorphismGroup],
function(g, G)
  local i;

  if IsOne(g) then return true;
  elif IsTrivial(G) then return false; fi;

  if CanEasilyTestSphericalTransitivity(G) and
    CanEasilyTestSphericalTransitivity(g) and
      not IsSphericallyTransitive(G) and
        IsSphericallyTransitive(g)
  then
    Info(InfoAutomGrp, 3, "g in G: false");
    Info(InfoAutomGrp, 3, "  not IsSphericallyTransitive(G) and IsSphericallyTransitive(g)");
    Info(InfoAutomGrp, 3, "  g = ", g);
    Info(InfoAutomGrp, 3, "  G = ", G);
    return false;
  fi;

  if not AbelImage(g) in AbelImage(G) then
    Info(InfoAutomGrp, 3, "g in G: false");
    Info(InfoAutomGrp, 3, "  not AbelImage(g) in AbelImage(G)");
    Info(InfoAutomGrp, 3, "  g = ", g);
    Info(InfoAutomGrp, 3, "  G = ", G);
    return false;
  fi;

  G := AG_ReducedForm(G);
  g := AG_ReducedForm(g);

  if FindElement(G, function(el) return el=g; end ,true, 8)<>fail then
    Info(InfoAutomGrp, 3, "g in G: FindElement returned true");
    return true;
  fi;

#TODO
  for i in [1..10] do
    if not PermOnLevel(g, i) in PermGroupOnLevel(G, i) then
      Info(InfoAutomGrp, 3, "g in G: false");
      Info(InfoAutomGrp, 3, "  not PermOnLevel(g, ",i,") in PermGroupOnLevel(G, ",i,")");
      Info(InfoAutomGrp, 3, "  g = ", g);
      Info(InfoAutomGrp, 3, "  G = ", G);
      return false;
    fi;
  od;

  TryNextMethod();
end);

InstallMethod(IsomorphismPermGroup, "for [IsTreeAutomorphismGroup and IsSelfSimilar]",
              [IsTreeAutomorphismGroup and IsSelfSimilar], SUM_FLAGS,
function(G)
  local H, F, gens_in_freegrp, pi, pi_bar, hom_function, inv_hom_function;
  H := PermGroupOnLevel(G, LevelOfFaithfulAction(G));
  return AG_GroupHomomorphismByImagesNC(G, H, GeneratorsOfGroup(G), GeneratorsOfGroup(H));
end);

InstallMethod(ContainsSphericallyTransitiveElement, "for [IsTreeAutomorphismGroup and IsSelfSimilar and IsActingOnBinaryTree]",
              [IsTreeAutomorphismGroup and IsSelfSimilar and IsActingOnBinaryTree],
function(G)
  local abels, indices, ab, trans;
  abels := List( GeneratorsOfGroup(G), AbelImage);
  if Length(abels)=0 then return false; fi;
  ab:=abels[1];
  trans:=One(ab)/(One(ab)+IndeterminateOfUnivariateRationalFunction(ab));
  for indices in Combinations([1..Length(abels)]) do
    if Sum(abels{indices})=trans then
      SetSphericallyTransitiveElement(G,Product(GeneratorsOfGroup(G){indices}));
      return true;
    fi;
  od;
  SetSphericallyTransitiveElement(G,fail);
  return false;
end);

InstallMethod(SphericallyTransitiveElement, "for [IsTreeAutomorphismGroup and IsSelfSimilar and IsActingOnBinaryTree]",
              [IsTreeAutomorphismGroup and IsSelfSimilar and IsActingOnBinaryTree],
function(G)
  ContainsSphericallyTransitiveElement(G);
  return SphericallyTransitiveElement(G);
end);

###############################################################################
##
#M  MarkovOperator(<G>, <lev>, <weights>)
##
InstallMethod(MarkovOperator, "for [IsTreeAutomorphismGroup, IsCyclotomic, IsList]",
              [IsTreeAutomorphismGroup, IsCyclotomic, IsList],
function(G, n, weights)
  local gens,R,indet;
  gens := ShallowCopy(GeneratorsOfGroup(G));
  if Length(weights)<>2*Length(gens) then Error("The number of weights must be twice as big as the number of generators"); fi;
  Append(gens, List(gens, x -> x^-1));
  if IsString(weights[1]) then
    R:=PolynomialRing(Integers,weights);
    indet:=IndeterminatesOfPolynomialRing(R);
    return Sum(List([1..Length(gens)], x -> indet[x]*PermOnLevelAsMatrix(gens[x], n)));
  else
    return Sum(List([1..Length(gens)], x -> weights[x]*PermOnLevelAsMatrix(gens[x], n)));
  fi;
end);

###############################################################################
##
#M  MarkovOperator(<G>, <lev>)
##
InstallMethod(MarkovOperator, "for [IsTreeAutomorphismGroup, IsCyclotomic]",
              [IsTreeAutomorphismGroup, IsCyclotomic],
function(G, n)
  return MarkovOperator(G,n,List([1..2*Length(GeneratorsOfGroup(G))],x->1/(2*Length(GeneratorsOfGroup(G)))));
end);

###############################################################################
##
#M  Section(<G>, <v>)
##
InstallMethod(Section, "for [IsTreeAutomorphismGroup, IsList]",
              [IsTreeAutomorphismGroup, IsList],
function(G, v)
  local gens, sec_gens, g;
  gens:=GeneratorsOfGroup(G);
  for g in gens do
    if v^g<>v then
      Error("Section([IsTreeAutomorphismGroup, IsList]): <G> does not fix <v>");
      return fail;
    fi;
  od;

  if not IsAutomGroup(G) then
    return Group(List(gens,x->Section(x,v)));
  else
    sec_gens:=[];
    for g in List(gens,x->Section(x,v)) do
      if not IsOne(g) and not g in sec_gens then
        Add(sec_gens, g);
      fi;
    od;
    if Length(sec_gens)>0 then
      return Group(sec_gens);
    else
      return Group([One(FamilyObj(Section(gens[1],v)))]);
    fi;
  fi;
end);

###############################################################################
##
#M  Section(<G>, <n>)
##
InstallMethod(Section, "for [IsTreeAutomorphismGroup, IsPosInt]",
              [IsTreeAutomorphismGroup, IsPosInt],
function(G, n)
  return Section(G,[n]);
end);

#E

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

SIMD_VMX.CellFixProblems := function(c)
    local s;
    s := UnionSet(Set(List(Collect(c, [@(0,[vload_2h_4x32f, vload_2l_4x32f]), @1.cond(e->ObjId(e)=var), @2]), i->i.args[1])),
         Set(List(Collect(c, [@(0,[vload1_8x16i, vload2_8x16i, vload4_8x16i, vloadu_8x16i]), @1.cond(e->ObjId(e)=var), @2, @(3)]), i->i.args[1])));

    if Length(s) > 0 then
        return decl(List(s), c);
    else
        return c;
    fi;
end;

ext: gi


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

# q:=HashLookup(bf, DFT(64))[1].ruletree;;
# s:= RCTDA(RC(SumsRuleTree(q)));;
# c:= _CodeSums(s, Y, X);;

# cc:=Copy(c);;
# cc:=Compile.pullDataAndDecls(cc);;
# cc:=Compile.fastScalarize(cc);;
# cc:=FlattenCode(UnrollCode(cc));;

# cc := CopyPropagate(cc);;

# IsLoop := c -> Same(ObjId(c),loop);

Declare(myUnrollLoop);

_myUnrollCode := function(c,decls)
     local op,d,len,tentry;
     op := ObjId(c);
     if Same(op,loop) then return myUnrollLoop(c,  decls);
     elif Same(op,chain) then return Concatenation(List(c.cmds, x->_myUnrollCode(x,decls)));
     elif Same(op,decl) then
	 Append(decls, c.vars);
	 for d in c.vars do
	    if IsArray(d.t) then
		d.value := V(List([1..d.t.size], i -> var.fresh_t("s", d.t.t))); fi;
	 od;
	 return _myUnrollCode(c.cmd,  decls);
     elif Same(op,data) then
	 #datas.(c.var.id) := c.value;
	 c.var.value := c.value;
	 return _myUnrollCode(c.cmd,  decls);
     elif Same(op,assign) then return [c];
     else Error("Can't handle object of type '", op, "'");
     fi;
end;

myUnrollLoop := function(c,  decls)
    local result, cmd, cmds, ind, val, i;
    cmds := _myUnrollCode(c.cmd,decls);
    ind := c.var;
    result := [];
    for i in c.range do
       val := V(i);
       Add(result, assign(ind, val));
       Append(result, Copy(cmds));
    od;
    return result;
end;

myUnrollCode := c -> _myUnrollCode(c,[]);


text:
function [index,distance]=near(x,x0,n);
% NEAR  finds the indices of x that are closest to the point x0.
%function [index,distance]=near(x,x0,[n]);
%     x is an array, x0 is a point, n is the number of closest points to get
%     (in order of increasing distance).  Distance is the abs(x-x0)
% rsignell@usgs.gov

if nargin==2
     n=1;
end
[distance,index]=sort(abs(x-x0));
distance=distance(1:n);
index=index(1:n);

end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function data=locsmooth(data,Fs,Tw,Ts)
%  Running line fit (using local linear regression) - 1d only, continuous
%  processes
%  Usage: data=locsmooth(data,Fs,Tw,Ts)
%  Inputs:
% Note that units of Fs, movinwin have to be consistent.
%  data  (single vector)
%  Fs    (sampling frequency) - optional. Default 1
%  Tw    (length of moving window) - optional.  Default. full length of data (global detrend)
%  Ts    (step size) - optional. Default Tw/2.
%
% Output:
% data   (locally smoothed data).
data=change_row_to_column(data);
N=size(data,1);
if nargin < 2; Fs=1; end;
if nargin < 3; Tw=N/Fs; end;
if nargin < 4; Ts=Tw/2; end;

n=round(Fs*Tw);
dn=round(Fs*Ts);
if ~isreal(data)
  yr=real(data);
  yi=imag(data);
  tmp=runline(yr,n,dn);
  yr=tmp';
  tmp=runline(yi,n,dn);
  yi=tmp';
  data=yr+i*yi;
else
  tmp=runline(data,n,dn);
  data=tmp';
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [y,fs]=readflac(filename,mode)
%READWAV  Read a .FLAC format sound file [Y,FS]=(FILENAME,MODE)
%
% Input Parameters:
%
%	FILENAME gives the name of the file (with optional .WAV extension) or alternatively
%                 can be the FIDX output from a previous call to READWAV
%	MODE		specifies the following (*=default):
%
%    Scaling: 's'    Auto scale to make data peak = +-1
%             'r'    Raw unscaled data (integer values)
%             'q'    Scaled to make 0dBm0 be unity mean square
%             'p' *	 Scaled to make +-1 equal full scale
%             'o'    Scale to bin centre rather than bin edge (e.g. 127 rather than 127.5 for 8 bit values)
%                     (can be combined with n+p,r,s modes)
%             'n'    Scale to negative peak rather than positive peak (e.g. 128.5 rather than 127.5 for 8 bit values)
%                     (can be combined with o+p,r,s modes)
%
% FLAC (Free lossless audio codec) is a compressed audio file format described here:
% http://flac.sourceforge.net/

%      Copyright (C) Mike Brookes 2008
%      Version: $Id: readflac.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2
    mode='p';
else
    mode = [mode(:).' 'p'];
end
dirt=voicebox('dir_temp');
[fnp,fnn,fne,fnv]=fileparts(filename);
filetemp=fullfile(dirt,[fnn '.wav' fnv]);
doscom=['"' voicebox('flac') '"' ' -d -f -o "' filetemp '" "' filename '"'];
%                     fprintf(1,'Executing: %s\n',doscom);
[doss,dosr]=dos(doscom); % run the program
if doss % test for errors
    error(sprintf('Error running DOS command: %s',doscom));
end
if exist(filetemp)~=2
    error(sprintf('No output file from: %s',doscom));
end
[y,fs]=readwav(filetemp,mode);
doscom=['del /f "' filetemp '"'];
if dos(doscom) % run the program
    error(sprintf('Error running DOS command: %s',doscom));
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
% pop_prepPipeline() - runs the early stage pipeline to reference and to
% detect bad channels
%
% Usage:
%   >>   [OUTEEG, com] = pop_prepPipeline(INEEG, params);
%
% Inputs:
%   INEEG   - input EEG dataset
%   params  - (optional) structure with parameters to override defaults
%
% Outputs:
%   OUTEEG  - output dataset
%
% See also:
%   prepPipeline, prepPipelineReport, EEGLAB

% Copyright (C) 2015  Kay Robbins with contributions from Nima
% Bigdely-Shamlo, Christian Kothe, Tim Mullen, and Cassidy Matousek
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [EEG, com] = pop_prepPipeline(EEG, params)
com = ''; % Return something if user presses the cancel button
okay = true;
if nargin < 1  %% display help if not enough arguments
    help pop_prepPipeline;
    return;
elseif nargin < 2
    params = struct();
end

%% Add path to prepPipeline subdirectories if not in the list
tmp = which('getPrepDefaults');
if isempty(tmp)
    myPath = fileparts(which('prepPipeline'));
    addpath(genpath(myPath));
end

%% Pop up window
if nargin < 2
    userData = getUserData();
    [params, okay] = MasterGUI([],[],userData, EEG);
end

%% Return if user pressed cancel
if (~okay)
	return;
end

%% Being the pipeline execution
userData = getUserData();
com = sprintf('%s = pop_prepPipeline(%s, %s);', inputname(1), ...
    struct2str(params));

reportMode = userData.report.reportMode.value;
consoleFID = userData.report.consoleFID.value;
publishOn = userData.report.publishOn.value;
summaryFilePath = userData.report.summaryFilePath.value;
sessionFilePath = userData.report.sessionFilePath.value;

if okay
    if strcmpi(reportMode, 'normal') || strcmpi(reportMode, 'skipReport')
        EEG = prepPipeline(EEG, params);
    end
    if (strcmpi(reportMode, 'normal') || ...
            strcmpi(reportMode, 'reportOnly')) %&& publishOn
        publishPrepReport(EEG, summaryFilePath, sessionFilePath, ...
            consoleFID, publishOn);
    end
    if strcmpi(reportMode, 'normal') || strcmpi(reportMode, 'skipReport')
        EEG = prepPostProcess(EEG, params);
    end
end

    function userData = getUserData()
        %% Gets the userData defaults and merges it with the parameters
        userData = struct('boundary', [], 'detrend', [], ...
            'lineNoise', [], 'reference', [], ...
            'report', [],  'postProcess', []);
        stepNames = fieldnames(userData);
        for k = 1:length(stepNames)
            defaults = getPrepDefaults(EEG, stepNames{k});
            [theseValues, errors] = checkStructureDefaults(params, ...
                defaults);
            if ~isempty(errors)
                error('pop_prepPipeline:BadParameters', ['|' ...
                    sprintf('%s|', errors{:})]);
            end
            userData.(stepNames{k}) = theseValues;
        end
    end  % getUserData

end % pop_prepPipeline

author: VisLab

repo: EEG-Clean-Tools


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setting GPU & Enviroment etc.
addpath('caffe/matlab');
caffe.reset_all();

caffe.set_mode_gpu();
gpu_id = 1;
caffe.set_device(gpu_id);
pp=0;

[foldername, videoname, filename, labelid] = textread('lst/UCF-101_testSplit1_25frames.lst','%s %s %d %d');

net_model = ['prototxt/spatial_cls.prototxt'];
net_weights = ['models/spatial_UCF101_list01.caffemodel'];
phase = 'test';
if ~exist(net_weights, 'file')
  error('Please download CaffeNet from Model Zoo before you run this demo');
end

net = caffe.Net(net_model, net_weights, phase);
im_data_ori = zeros(340,256, 3,25,'single');
mean_data = caffe.io.read_mean('mean/VGG_mean.binaryproto');

input_data = zeros(224,224,3,25);
scores_total=zeros(101,size(foldername,1)/25);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Processing videos
for i=1:25:size(filename)
    ii=1;
    t_io = 0;
    t1 = clock;
    for j=0:24
	for k=filename(i+j):filename(i+j)
	  i_str=sprintf('%04d',k);
      t2 = clock;
      im = caffe.io.load_image(['ucf101_image/',char(foldername{i+j}),'/',char(videoname{i+j}),'/flow_i_',i_str,'.jpg']);
      t3 = clock;
      t_io = etime(t3,t2) + t_io;

	  im_data_ori(:,:,:,j+1) = im(:,:,:);
%      im_data_mirror(:,:,:,j+1) = im(end:-1:1,:,:);
	end
    end

	input_data = im_data_ori(60:60+223,16:16+223,:,:);
	input_data = bsxfun(@minus,input_data,mean_data);

    net.blobs('data').set_data(input_data);
    net.forward_prefilled();
    scores = net.blobs('prob').get_data();
    scores_total(:,(i-1)/25+1) = mean(scores, 2);
    t4 = clock;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    total_time = etime(t4,t1);
    io_time = t_io;
    computational_time = etime(t4,t1)-t_io;
    fprintf(['Time (w/o I/O) of processing 25 rgb frames: ',num2str(computational_time),'\n'])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%show results:
    [~, maxlabel] = max(mean(scores,2));
    if (maxlabel==labelid(i)+1)
        pp=pp+1;
    end
    fprintf('%d/%d %d %d %f\n',pp,(i-1)/25+1,maxlabel, labelid(i)+1, pp/((i-1)/25+1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cleaning and etc.
save('UCF101_spatial_list01.mat','scores_total');
caffe.reset_all();
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: zbwglory

repo: MV-release

save_path: github-repos/MATLAB/zbwglory-MV-release


text:
classdef ROIPooling < dagnn.Layer
  % DAGNN.ROIPOOLING  Region of interest pooling layer

  % Copyright (C) 2016 Hakan Bilen.
  % All rights reserved.
  %
  % This file is part of the VLFeat library and is made available under
  % the terms of the BSD license (see the COPYING file).

  properties
    method = 'max'
    subdivisions = [6 6]
    transform = 1
    flatten = false
  end

  methods
    function outputs = forward(obj, inputs, params)
      numROIs = numel(inputs{2}) / 5 ;
      outputs{1} = vl_nnroipool(...
        inputs{1}, inputs{2}, ...
        'subdivisions', obj.subdivisions, ...
        'transform', obj.transform, ...
        'method', obj.method) ;
      if obj.flatten
        outputs{1} = reshape(outputs{1},1,1,[],numROIs) ;
      end
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      numROIs = numel(inputs{2}) / 5 ;
      if obj.flatten
        % unflatten
        derOutputs{1} = reshape(...
          derOutputs{1},obj.subdivisions(1),obj.subdivisions(2),[],numROIs) ;
      end
      derInputs{1} = vl_nnroipool(...
        inputs{1}, inputs{2}, derOutputs{1}, ...
        'subdivisions', obj.subdivisions, ...
        'transform', obj.transform, ...
        'method', obj.method) ;
      derInputs{2} = [];
      derParams = {} ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      if isempty(inputSizes{1})
        n = 0 ;
      else
        n = prod(inputSizes{2})/5 ;
      end
      outputSizes{1} = [obj.subdivisions, inputSizes{1}(3), n] ;
    end

    function obj = ROIPooling(varargin)
      obj.load(varargin) ;
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function handles = togglePlotSilhouette(hObject, handles, state)

set(handles.btnPlotSilhouette, 'Enable', state);

if isequal(handles.foundExportFig, 1)
    set(handles.btnPlotSilhouetteSave, 'Enable', state);
end

guidata(hObject, handles);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function xlsPasteTo(filename,sheetname,width, height,varargin)

%Paste current figure to selected Excel sheet and cell
%
%
% xlsPasteTo(filename,sheetname,width, height,range)
%Example:
%xlsPasteTo('File.xls','Sheet1',200, 200,'A1')
% this will paset into A1 at Sheet1 at File.xls the current figure with
% width and height of 200
%
% tal.shir@hotmail.com

options = varargin;

    range = varargin{1};

[fpath,file,ext] = fileparts(char(filename));
if isempty(fpath)
    fpath = pwd;
end
Excel = actxserver('Excel.Application');
set(Excel,'Visible',0);
Workbook = invoke(Excel.Workbooks, 'open', [fpath filesep file ext]);

sheet = get(Excel.Worksheets, 'Item',sheetname);
invoke(sheet,'Activate');

    ExAct = Excel.Activesheet;

  ExActRange = get(ExAct,'Range',range);
    ExActRange.Select;
    pos=get(gcf,'Position');
    set(gcf,'Position',[ pos(1:2) width height])
    print -dmeta

invoke(Excel.Selection,'PasteSpecial');

invoke(Workbook, 'Save');
invoke(Excel, 'Quit');
delete(Excel);

author: Sable

repo: mcbench-benchmarks


text:
%% Copyright (C) 2016-2017, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym ezsurf (@var{z})
%% @defmethodx @@sym ezsurf (@var{f1}, @var{f2}, @var{f3})
%% @defmethodx @@sym ezsurf (@dots{}, @var{dom})
%% @defmethodx @@sym ezsurf (@dots{}, @var{N})
%% Simple 3D surface plots of symbolic expressions.
%%
%% Example 3D surface plot:
%% @example
%% @group
%% syms x y
%% z = sin(2*x)*sin(y)
%%   @result{} z = (sym) sin(2⋅x)⋅sin(y)
%% ezsurf(z)                                    % doctest: +SKIP
%% @end group
%% @end example
%%
%% Example parametric surface plot of a Möbius strip:
%% @example
%% @group
%% syms u v
%% x = (1+v*cos(u/2))*cos(u)
%%   @result{} x = (sym)
%%       ⎛     ⎛u⎞    ⎞
%%       ⎜v⋅cos⎜─⎟ + 1⎟⋅cos(u)
%%       ⎝     ⎝2⎠    ⎠
%% y = (1+v*cos(u/2))*sin(u);
%% z = v*sin(u/2);
%%
%% ezsurf(x, y, z, [0 2*pi -0.5 0.5], 32)       % doctest: +SKIP
%% axis equal
%% @end group
%% @end example
%%
%% See help for the (non-symbolic) @code{ezsurf}, which this
%% routine calls after trying to convert sym inputs to
%% anonymous functions.
%%
%% @seealso{ezsurf, @@sym/ezmesh, @@sym/ezplot, @@sym/function_handle}
%% @end defmethod

function varargout = ezsurf(varargin)

  % first input is handle, shift
  if (ishandle(varargin{1}))
    firstpotsym = 2;
  else
    firstpotsym = 1;
  end

  maxnumsym = 3;
  firstsym = [];

  for i = firstpotsym:nargin
    if (isa(varargin{i}, 'sym'))
      if (i < firstpotsym + maxnumsym)
        % one of the fcns to plot, covert to handle fcn

        % Each is function of one var, and its the same var for all
        thissym = symvar(varargin{i});
        assert(length(thissym) <= 2, ...
          'ezsurf: parameterized: functions should have at most two inputs');
        if (isempty(thissym))
          % a number, create a constant function in a dummy variable
          % (0*t works around some Octave oddity on 3.8 and hg Dec 2014)
          thisf = inline(sprintf('%g + 0*t', double(varargin{i})), 't');
          %thisf = @(t) 0*t + double(varargin{i});  % no
        else
          % check variables match (sanity check)
          if (isempty(firstsym))
            firstsym = thissym;
          else
            assert(all(logical(thissym == firstsym)), ...
              'ezsurf: all functions must be in terms of the same variables');
          end
          thisf = function_handle(varargin{i});
        end

        varargin{i} = thisf;

      else
        % plot ranges, etc, convert syms to doubles
        varargin{i} = double(varargin{i});
      end
    end
  end

  h = ezsurf(varargin{:});

  if (nargout)
    varargout{1} = h;
  end

end

%!error <all functions must be in terms of the same variables>
%! syms u v t
%! ezsurf(u*v, 2*u*v, 3*v*t)

%!error <functions should have at most two inputs>
%! syms u v t
%! ezsurf(u*v, 2*u*v, u*v*t)

repo: octsympy


text:
function [ xc, rank, seed ] = comp_random_grlex ( kc, rank1, rank2, seed )

%*****************************************************************************80
%
%% COMP_RANDOM_GRLEX randomly selects a composition in a given rank range.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer KC, the number of parts of the composition.
%    1 <= KC.
%
%    Input, integer RANK1, RANK2, the lowest and highest ranks to consider.
%    1 <= RANK1 <= RANK2.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
%    Output, integer XC(KC), a randomly chosen composition with rank
%    between RANK1 and RANK2.
%
%    Output, integer RANK, the rank of the composition.
%

%
%  Ensure that 1 <= KC.
%
  if ( kc < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'COMP_RANDOM_GRLEX - Fatal error!\n' );
    fprintf ( 1, '  KC < 1\n' );
    error ( 'COMP_RANDOM_GRLEX - Fatal error!' );
  end
%
%  Ensure that 1 <= RANK1.
%
  if ( rank1 < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'COMP_RANDOM_GRLEX - Fatal error!\n' );
    fprintf ( 1, '  RANK1 < 1\n' );
    error ( 'COMP_RANDOM_GRLEX - Fatal error!' );
  end
%
%  Ensure that RANK1 <= RANK2.
%
  if ( rank2 < rank1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'COMP_RANDOM_GRLEX - Fatal error!\n' );
    fprintf ( 1, '  RANK2 < RANK1\n' );
    error ( 'COMP_RANDOM_GRLEX - Fatal error!' );
  end
%
%  Select the rank.
%
  [ rank, seed ] = i4_uniform_ab ( rank1, rank2, seed );
%
%  Recover the composition.
%
  xc = comp_unrank_grlex ( kc, rank );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% SHOW - Displays an image with the right size and colors and with a title.
%
% Usage:
%         h = show(im)
%         h = show(im, figNo)
%         h = show(im, title)
%         h = show(im, figNo, title)
%
% Arguments:  im    - Either a 2 or 3D array of pixel values or the name
%                     of an image file;
%             figNo - Optional figure number to display image in. If
%                     figNo is 0 the current figure or subplot is
%                     assumed.
%             title - Optional string specifying figure title
%
% Returns:    h     - Handle to the figure.  This allows you to set
%                     additional figure attributes if desired.
%
% The function displays the image, automatically setting the colour map to
% grey if it is a 2D image, or leaving it as colour otherwise, and setting
% the axes to be 'equal'.  The image is also displayed as 'TrueSize', that
% is, pixels on the screen match pixels in the image (if it is possible
% to fit it on the screen, otherwise MATLAB rescales it to fit).
%
% Unless you are doing a subplot (figNo==0) the window is sized to match
% the image, leaving no border, and hence saving desktop real estate.
%
% If figNo is omitted a new figure window is created for the image.  If
% figNo is supplied, and the figure exists, the existing window is reused to
% display the image, otherwise a new window is created. If figNo is 0 the
% current figure or subplot is assumed.
%
% See also: SHOWSURF

% Copyright (c) 2000-2009 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% October   2000  Original version
% March     2003  Mods to alow figure name in window bar and allow for subplots.
% April     2007  Proper recording and restoring of MATLAB warning state.
% September 2008  Octave compatible
% May       2009  Reworked argument handling logic for extra flexibility

function h = show(im, param2, param3)

    Octave = exist('OCTAVE_VERSION') ~= 0;  % Are we running under Octave?

    if ~Octave
	s = warning('query','all'); % Record existing warning state.
	warning('off');             % Turn off warnings that might arise if image
                                    % has to be rescaled to fit on screen
    end

    % Check case where im is an image filename rather than image data
    if ~isnumeric(im) & ~islogical(im)
        Title = im;            % Default title is file name
        im = imread(im);
    else
        Title = inputname(1);  % Default title is variable name of image data
    end

    figNo = -1;                % Default value indicating create new figure

    % If two arguments check type of 2nd argument to see if it is the title or
    % figure number that has been supplied
    if nargin == 2
        if strcmp(class(param2),'char')
            Title = param2;
        elseif isnumeric(param2) && length(param2) == 1
            figNo = param2;
        else
            error('2nd argument must be a figure number or title');
        end
    elseif nargin == 3
        figNo = param2;
        Title = param3;

        if ~strcmp(class(Title),'char')
            error('Title must be a string');
        end

        if ~isnumeric(param2) || length(param2) ~= 1
            error('Figure number must be an integer');
        end
    end

    if figNo > 0           % We have a valid figure number
        figure(figNo);     % Reuse or create a figure window with this number
        if ~Octave
            subplot('position',[0 0 1 1]); % Use the whole window
        end
    elseif figNo == -1
        figNo = figure;        % Create new figure window
        if ~Octave
            subplot('position',[0 0 1 1]); % Use the whole window
        end
    end

    if ndims(im) == 2          % Display as greyscale
	imagesc(im);
	colormap('gray');
    else
	imshow(im);            % Display as RGB
    end

    if figNo == 0                 % Assume we are trying to do a subplot
        figNo = gcf;              % Get the current figure number
	axis('image'); axis('off');
	title(Title);             % Use a title rather than rename the figure
    else
	axis('image'); axis('off');
	set(figNo,'name', ['  ' Title])
	if ~Octave
	    truesize(figNo);
	end
    end

    if nargout == 1
       h = figNo;
    end

    if ~Octave
	warning(s);  % Restore warnings
    end

author: jianxiongxiao

repo: ProfXkit


text:
function [resdQ, rhsQ] = volumeAveragedConcentrationFlux(dQ,Q,jflux,T,param)
%   volumeAveragedConcentrationFlux is used to implement the three parameters reduced model for solid phase diffusion.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% Diffusion coefficients for the solid phase
[Dps_eff, Dns_eff] = param.SolidDiffusionCoefficientsFunction(T,param);

rhsQ_p      = -30*Dps_eff./param.Rp_p^2.*Q(1:param.Np)    -45./(2*param.Rp_p^2).*jflux(1:param.Np);

resdQ_p     = dQ(1:param.Np)        -rhsQ_p;

rhsQ_n      = -30*Dns_eff./param.Rp_n^2.*Q(param.Np+1:end)-45./(2*param.Rp_n^2).*jflux(param.Np+1:end);

resdQ_n     = dQ(param.Np+1:end) - rhsQ_n   ;

rhsQ        = [rhsQ_p;rhsQ_n];
resdQ       = [resdQ_p;resdQ_n];

% This model has been taken from the paper, "Efficient Macro-Micro Scale Coupled
% Modeling of Batteries" -  Subramanian,Diwakar,Tapriyal - 2005 JES

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function e = end(T,i,n)
% overloaded end function

if n==1
  e = prod(size(T)); %#ok<PSIZE>
else
  e = size(T,i);
end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [ obj ] = vb_invert( obj )
% Run VB update equations until convergence.
%
% This is a protected method of the tapas_Huge class. It cannot be called
% from outside the class.
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch)
% Copyright (C) 2019 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This software is provided "as is", without warranty of any kind, express
% or implied, including, but not limited to the warranties of
% merchantability, fitness for a particular purpose and non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

%% initialization
obj = obj.vb_init( );
bUpdateKmeans       = false;
bUpdateAssignments  = false;
bUpdateClusters     = false;
bVerbose            = obj.options.nvp.verbose;
nIt = obj.options.nvp.numberofiterations;
if isempty(nIt)
    nIt = 999;
end
obj.trace.nfe = repmat(obj.trace.nfe, nIt, 1);

%% ======= MAIN LOOP =======
for iIt = 1:nIt
% schedule posterior updates

    if bUpdateKmeans
        % initialize cluster parameters with kmeans
        obj = kmeans_init( obj );
        bUpdateKmeans = false;
    end

    if bUpdateAssignments
        % update alpha, q_nk
        obj = update_assigments( obj );
    end

    if bUpdateClusters
        if obj.options.confVar
            % update m_beta, S_beta
            obj = update_confounds( obj );
        end
        % update m_k, tau_k, nu_k, S_k
        obj = update_clusters( obj ); %%% TODO add support for cluster-specific prior mean and scale
    end

    % update mu_n, Sigma_n, b_nr
    obj = update_dcm( obj, bUpdateClusters );

    % update negative free energy
    obj.trace.nfe(iIt) = obj.vb_nfe( );
    dNfe = obj.trace.nfe(iIt) - obj.posterior.nfe;
    obj.posterior.nfe = obj.trace.nfe(iIt);
    if bVerbose
        if isnan(dNfe)
            fprintf('iteration %3u.\n', iIt);
        else
            fprintf('iteration %3u, dF: %.2E.\n', iIt, dNfe);
        end
    end

    % check stopping conditions
    if ~bUpdateClusters
        % start updating clustering parameters and confound coefficients
        bUpdateClusters = (dNfe < obj.options.convergence.dDcm && dNfe >= 0);% || iIt > 32;
        if bUpdateClusters
            obj.trace.convergence(1) = iIt;
        end

    elseif ~bUpdateAssignments
        % start updating assignment labels
        bUpdateAssignments = dNfe < obj.options.convergence.dClusters;
        if bUpdateAssignments
            bUpdateKmeans = obj.options.convergence.bUseKmeans;
            obj.trace.convergence(2) = iIt;
            % force next iteration to happen
            obj.posterior.nfe = NaN;
        end

    elseif dNfe < obj.options.convergence.dVb
        % convergence of VB scheme
        obj.trace.convergence(3) = iIt;
        obj.trace.nfe(iIt+1:end) = [];
        obj.trace.epsilon = obj.aux.epsilon;
        if bVerbose
            fprintf('VB converged after %u iterations. ', iIt)
        end
        break;

    end
end

if bVerbose
    if iIt == nIt
        fprintf('VB reached the maximum number of iterations. ')
    end
    fprintf('Negative free energy: %.2E\n', obj.posterior.nfe);
end

%% post-processing
% estimate confound-free DCM parameters
if obj.options.confVar
    mu_beta = obj.aux.mu_beta;
    Sigma_beta = obj.aux.Sigma_beta;
    if obj.options.confVar > 1
        mu_beta = sum(bsxfun(@times, mu_beta, permute(obj.posterior.q_nk, ...
            [1 3 2])), 3);
        Sigma_beta = sum(bsxfun(@times, Sigma_beta, permute( ...
            obj.posterior.q_nk, [3 4 1 2])), 4);
    end
    obj.posterior.mu_r = obj.posterior.mu_n(:, 1:obj.idx.P_c) - mu_beta;
    obj.posterior.Sigma_r = obj.posterior.Sigma_n(1:obj.idx.P_c, 1:obj.idx.P_c, :) ...
        + Sigma_beta;
else
    obj.posterior.mu_r = [];
    obj.posterior.Sigma_r = [];
end

% normalized variance of residual (by subject and region)
obj.posterior.nrv = zeros(obj.N, obj.R);
for n = 1:obj.N
    varRes = var(obj.aux.epsilon{n});
    varSig = var(obj.data(n).bold);
    obj.posterior.nrv(n,:) = varRes./varSig;
end

end

%% UPDATE: posterior over cluster parameters
function [ obj ] = update_clusters( obj )

% update parameters of inverse Wishart
obj.posterior.tau = obj.prior.tau_0 + obj.aux.q_k;
obj.posterior.nu  = obj.prior.nu_0  + obj.aux.q_k;

% auxiliary variables
Sigma_n_c = obj.posterior.Sigma_n(1:obj.idx.P_c,1:obj.idx.P_c,:);
mu_n_c = obj.posterior.mu_n(:,1:obj.idx.P_c);
if obj.options.confVar == 1
    Sigma_n_c = Sigma_n_c + obj.aux.Sigma_beta;
    mu_n_c = mu_n_c - obj.aux.mu_beta;
end

for k = 1:obj.K

    if obj.options.confVar > 1
        Sigma_n_c = obj.posterior.Sigma_n(1:obj.idx.P_c,1:obj.idx.P_c,:) + ...
            obj.aux.Sigma_beta(:,:,:,k);
        mu_n_c = obj.posterior.mu_n(:,1:obj.idx.P_c) - ...
            obj.aux.mu_beta(:,:,k);
    end
    % Eq (17)
    mu_k_c = obj.posterior.q_nk(:,k)'*mu_n_c;
    Sigma_k_c = sum(bsxfun(@times, Sigma_n_c, ...
        permute(obj.posterior.q_nk(:,k),[3,2,1])), 3);
    % Eq (16)
    % posterior cluster mean
    obj.posterior.m(k,:) = (mu_k_c + obj.prior.tau_0*obj.prior.m_0)/...
        (obj.aux.q_k(k) + obj.prior.tau_0);
    % posterior cluster covariance
    mu_k_c = mu_k_c./obj.aux.q_k(k);
    dmu = bsxfun(@minus, mu_n_c, mu_k_c);
    obj.posterior.S(:,:,k) = ...
        Sigma_k_c + ...
        dmu'*bsxfun(@times, dmu, obj.posterior.q_nk(:,k)) + ...
        obj.aux.q_k(k)*obj.prior.tau_0/...
            (obj.aux.q_k(k) + obj.prior.tau_0)*...
            (mu_k_c - obj.prior.m_0)'*(mu_k_c - obj.prior.m_0) + ...
        obj.prior.S_0;
    % update auxiliary variables
    obj.aux.ldS(k) = tapas_huge_logdet(obj.posterior.S(:,:,k));
    obj.aux.nu_inv_S(:,:,k) = obj.posterior.nu(k).*inv(obj.posterior.S(:,:,k));

end
end

%% UPDATE: posterior over cluster assignments
function [ obj ] = update_assigments( obj )

Sigma_n_c = obj.posterior.Sigma_n(1:obj.idx.P_c,1:obj.idx.P_c,:);
mu_n_c = obj.posterior.mu_n(:,1:obj.idx.P_c);
if obj.options.confVar == 1
    Sigma_n_c = Sigma_n_c + obj.aux.Sigma_beta;
    mu_n_c = mu_n_c - obj.aux.mu_beta;
end
psiNu = sum(psi(0, bsxfun(@minus, obj.posterior.nu, 0:obj.idx.P_c-1)/2), 2);
psiAlpha = psi(0, obj.posterior.alpha);
log_q_nk = zeros(obj.N, obj.K);
% Eq (18)
for k = 1:obj.K

    if obj.options.confVar > 1
        Sigma_n_c = obj.posterior.Sigma_n(1:obj.idx.P_c,1:obj.idx.P_c,:) + ...
            obj.aux.Sigma_beta(:,:,:,k);
        mu_n_c = obj.posterior.mu_n(:,1:obj.idx.P_c) - ...
            obj.aux.mu_beta(:,:,k);
    end
    dm = bsxfun(@minus,mu_n_c,obj.posterior.m(k,:));
    log_q_nk(:,k) = -.5*sum((dm*obj.aux.nu_inv_S(:,:,k)).*dm, 2) ...
        -.5*squeeze(sum(sum(bsxfun(@times, Sigma_n_c, ...
            obj.aux.nu_inv_S(:,:,k)), 1), 2));

end
log_q_nk = bsxfun(@plus,log_q_nk, (-.5*obj.aux.ldS +.5*psiNu ...
    -.5*obj.idx.P_c./obj.posterior.tau + psiAlpha)');
log_q_nk = exp(bsxfun(@minus, log_q_nk, max(log_q_nk, [], 2)));
obj.posterior.q_nk = bsxfun(@rdivide, log_q_nk, sum(log_q_nk, 2));

% Eq (17)
obj.aux.q_k = sum(obj.posterior.q_nk,1)' + realmin;

% Eq (15)
obj.posterior.alpha = obj.prior.alpha_0 + obj.aux.q_k;

end

%% UPDATE: posterior over confound coefficients
function [ obj ] = update_confounds( obj )

m_k = permute(obj.posterior.m, [3,2,1]);
mu_n_c = obj.posterior.mu_n(:, 1:obj.idx.P_c);

for p = 1:obj.idx.P_c

    % Eq (XX)
    tmp = obj.aux.nu_inv_S(p, p, :);
    tmp = bsxfun(@times, obj.posterior.q_nk, tmp(:)');
    if obj.options.confVar == 1
        tmp = sum(tmp, 2);
    end
    % posterior precision over beta
    obj.posterior.S_beta(:, :, p, :) = bsxfun(@plus, obj.aux.Pi_beta_0, ...
        sum(bsxfun(@times, obj.aux.x_n_2, reshape(tmp, 1, 1, obj.N, [])), 3));

    % Eq (XX+1)
    obj.aux.mu_beta(:,p,:) = 0;
    dmu = bsxfun(@minus, mu_n_c, obj.aux.mu_beta);
    dmu = bsxfun(@minus, dmu, m_k);
    tmp = permute(obj.posterior.q_nk, [1 3 2]).*...
        sum(bsxfun(@times, obj.aux.nu_inv_S(p, :, :), dmu), 2);
    if obj.options.confVar == 1
        tmp = sum(tmp, 3);
    end
    tmp = permute(sum(bsxfun(@times, tmp, obj.aux.x_n), 1), [2 1 3]);
    for k = 1:size(obj.posterior.m_beta, 3)
        obj.posterior.m_beta(:,p,k) = obj.posterior.S_beta(:,:,p,k)\...
            (obj.aux.Pi_m_beta_0 + tmp(:,1,k));
        % posterior covariance over beta
        obj.posterior.S_beta(:,:,p,k) = inv(obj.posterior.S_beta(:,:,p,k));
        obj.aux.ldS_beta(p,k) = tapas_huge_logdet( ...
            obj.posterior.S_beta(:,:,p,k));
    end
    % update auxiliary variables
    obj.aux.mu_beta(:,p,:) = sum(bsxfun(@times, obj.aux.x_n, ...
        permute(obj.posterior.m_beta(:,p,:), [2 1 3])), 2);
    obj.aux.Sigma_beta(p,p,:,:) = sum(sum(bsxfun(@times, obj.aux.x_n_2, ...
        obj.posterior.S_beta(:,:,p,:)), 1), 2);

end

end

%% UPDATE: posterior over DCM and noise parameters
function [ obj ] = update_dcm( obj, bCheckNfe )

% cache DCM parameters
tmpMu       = obj.posterior.mu_n;
tmpSigma    = obj.posterior.Sigma_n;
tmpLdSigma  = obj.aux.ldSigma;
tmpEpsilon  = obj.aux.epsilon;
tmpG        = obj.aux.G;
tmpB        = obj.posterior.b;
tmpBp       = obj.aux.b_prime;
tmpLb       = obj.aux.lambda_bar;
tmpNfe      = obj.vb_nfe( );

% subject-wise update DCM and noise parameters
for n = 1:obj.N
    varEps = sum(var(obj.aux.epsilon{n}));
    % Eq (20)
    m_k = obj.posterior.m;
    if obj.options.confVar
        tmp = permute(obj.aux.mu_beta(n,:,:), [3 2 1]);
        m_k = bsxfun(@plus, m_k, tmp);
    end
    lambda_prime_c = zeros(obj.idx.P_c);
    mu_prime_c = zeros(1,obj.idx.P_c);
    for k = 1:obj.K
        tmps = obj.posterior.q_nk(n,k)*obj.aux.nu_inv_S(:,:,k);
        lambda_prime_c = lambda_prime_c + tmps;
        mu_prime_c = mu_prime_c + m_k(k,:)*tmps;
    end

    % Eq (19)
    % posterior covariance
    Lambda_bar = repmat(obj.aux.lambda_bar(n,:), obj.aux.q_r(n), 1);
    G_Lambda = bsxfun(@times, obj.aux.G{n}, Lambda_bar(:))';
    Pi_n = G_Lambda*obj.aux.G{n};
    Pi_n(1:obj.idx.P_c, 1:obj.idx.P_c) = lambda_prime_c + ...
        Pi_n(1:obj.idx.P_c, 1:obj.idx.P_c);
    Pi_n(obj.idx.P_c+1:end, obj.idx.P_c+1:end) = obj.aux.Pi_h + ...
        Pi_n(obj.idx.P_c+1:end, obj.idx.P_c+1:end);
    obj.posterior.Sigma_n(:,:,n) = inv(Pi_n);
    obj.aux.ldSigma(n) = - tapas_huge_logdet(Pi_n);
    % posterior mean
    tmp = obj.aux.epsilon{n}(:) + obj.aux.G{n}*obj.posterior.mu_n(n,:)';
    % regularization factor
    tmp = tmp*(.9.^mod(obj.trace.nRetract(n),5));
    obj.posterior.mu_n(n,:) = ...
        ((G_Lambda*tmp)' + [mu_prime_c, obj.aux.mu_prime_h])/Pi_n;
    % recalculate jacobian and residual
    obj = obj.options.fncBold(obj, n);

    % Eq (22)
    tmp = sum((obj.aux.G{n}*obj.posterior.Sigma_n(:,:,n)).*obj.aux.G{n},2);
    obj.aux.b_prime(n,:) = sum(obj.aux.epsilon{n}.^2 + ...
        reshape(tmp, obj.aux.q_r(n), obj.R), 1);
    % Eq (21)
    obj.posterior.b(n,:) = obj.prior.b_0 + obj.aux.b_prime(n,:)/2;
    % Eq (23)
    obj.aux.lambda_bar(n,:) = obj.posterior.a(n,:)./obj.posterior.b(n,:);

    newNfe = obj.vb_nfe( );
    % check negative free energy (during first few iterations, accept
    % decrease in F if fit improves)
    bUnstable = any(isnan(obj.aux.G{n}(:))) || any(isinf(obj.aux.G{n}(:)));
    bOverride = bCheckNfe || varEps < sum(var(obj.aux.epsilon{n}));
    if bUnstable || (tmpNfe > newNfe && bOverride)
        % retract update if new parameters unstable or if F decreased
        obj.posterior.mu_n(n,:)      = tmpMu(n,:);
        obj.posterior.Sigma_n(:,:,n) = tmpSigma(:,:,n);
        obj.aux.ldSigma(n)           = tmpLdSigma(n);
        obj.aux.epsilon{n}           = tmpEpsilon{n};
        obj.aux.G{n}                 = tmpG{n};
        obj.posterior.b(n,:)         = tmpB(n,:);
        obj.aux.b_prime(n,:)         = tmpBp(n,:);
        obj.aux.lambda_bar(n,:)      = tmpLb(n,:);
        obj.trace.nRetract(n)        = obj.trace.nRetract(n) + 1;
    else
        % otherwise accept update
        tmpNfe = newNfe;
        obj.trace.nDcmUpdate(n) = obj.trace.nDcmUpdate(n) + 1;
        obj.trace.nRetract(n)   = 0;
    end

end

end

%% initialize cluster parameters with K-means
function [ obj ] = kmeans_init( obj )
% kmeans
mu_n_c = obj.posterior.mu_n(:,1:obj.idx.P_c);
if obj.options.confVar
    mu_n_c = mu_n_c - obj.aux.mu_beta(:,:,1);
end
idx = kmeans(mu_n_c, obj.K, 'Replicates', obj.const.nKmeans);

% initialize additional clusters
obj.posterior.q_nk = zeros(obj.N,obj.K);
obj.posterior.q_nk(sub2ind([obj.N,obj.K], (1:obj.N)', idx)) = 1;
obj.aux.q_k = sum(obj.posterior.q_nk, 1)' + realmin;
obj.posterior.alpha = obj.prior.alpha_0 + obj.aux.q_k;
if obj.options.confVar > 1
    obj.posterior.m_beta = repmat(obj.posterior.m_beta(:,:,1), [1 1 obj.K]);
    obj.posterior.S_beta = repmat(obj.posterior.S_beta(:,:,:,1), ...
        [1 1 1 obj.K]);
    obj.aux.Sigma_beta = repmat(obj.aux.Sigma_beta(:,:,:,1), [1 1 1 obj.K]);
    obj.aux.ldS_beta = repmat(obj.aux.ldS_beta(:,1), 1, obj.K);
    obj.aux.mu_beta = repmat(obj.aux.mu_beta(:,:,1), [1 1 obj.K]);
end
obj = update_clusters( obj );

% debug info
obj.trace.kmeans.labels = idx;
obj.trace.kmeans.inputs = mu_n_c;

end

author: translationalneuromodeling

repo: tapas


text:
function MRSIStruct = setDimension(MRSIStruct, dimLabel, value)
    switch(lower(dimLabel))
        case 'x'
            MRSIStruct.dims.x = value;
        case 'y'
            MRSIStruct.dims.y = value;
        case 'z'
            MRSIStruct.dims.z = value;
        case 'coils'
            MRSIStruct.dims.coils = value;
        case 'averages'
            MRSIStruct.dims.averages = value;
        case 'kx'
            MRSIStruct.dims.kx = value;
        case 'ky'
            MRSIStruct.dims.ky = value;
        case 'kz'
            MRSIStruct.dims.kz = value;
        case 'subspecs'
            MRSIStruct.dims.subspecs = value;
        case 'extras'
            MRSIStruct.dims.extras = value;
        case 't'
            MRSIStruct.dims.t = value;
        otherwise
            error('Dimension not found %s', dimLabel);
    end
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function poly_fprintf ( file_name, node_num, segment, edge_num, ...
  edge_nodes, hole_num, hole_point )

%*****************************************************************************80
%
%% POLY_WRITE writes data to a POLY file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the name of the file to create.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real SEGMENT(2,NODE_NUM), the nodes.
%
%    Input, integer EDGE_NUM, the number of edges.
%
%    Input, integer EDGE_NODES(2,EDGE_NUM), the nodes that form each edge.
%
%    Input, integer HOLE_NUM, the number of holes in the mesh.
%
%    Input, real HOLE_POINT(2,HOLE_NUM), a point in each hole.
%
  region_num = 0;

  poly_unit = fopen ( file_name, 'w' );

  if ( poly_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'POLY_WRITE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the output POLY file.\n' );
    error ( 'POLY_WRITE - Fatal error!' );
  end

  fprintf ( poly_unit, '#  %s\n', file_name );
  fprintf ( poly_unit, '#  Created by poly_write.m\n' );
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '#  Vertex  Dimension  Attribute  Marker\n' );
  fprintf ( poly_unit, '#  Count              Count      0/1\n' );
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '  %d  2  0  0\n', node_num );
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '#  Vertex  X  Y  Attributes  Marker\n' );
  fprintf ( poly_unit, '#  Index\n' );
  fprintf ( poly_unit, '#\n' );
  for node = 1 : node_num
    fprintf ( poly_unit, '  %8d  %10f  %10f\n', node, segment(1:2,node) );
  end
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '#  Segment  Marker\n' );
  fprintf ( poly_unit, '#  Count    0/1\n' );
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '  %d  0\n', edge_num );
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '#  Segment_index  Node1  Node2  Marker\n' );
  fprintf ( poly_unit, '#\n' );
  for edge = 1 : edge_num
    fprintf ( poly_unit, '  %8d  %8d  %8d  0\n', edge, edge_nodes(1:2,edge) );
  end
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '#  Hole\n' );
  fprintf ( poly_unit, '#  Count\n' );
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '  %d\n', hole_num );
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '#  Hole_index  X  Y\n' );
  fprintf ( poly_unit, '#\n' );
  for hole = 1 : hole_num
    fprintf ( poly_unit, '  %8d  %10f  %10f\n', hole, hole_point(1:2,hole) );
  end
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '#  Region\n' );
  fprintf ( poly_unit, '#  Count\n' );
  fprintf ( poly_unit, '#\n' );
  fprintf ( poly_unit, '  %d\n', region_num );

  fclose ( poly_unit );

  return
end


text:
classdef VademecumLoader < handle

    properties (SetAccess = private, GetAccess = public)
        dataBase
    end

    properties (Access = private)
        index
        path
        iIndex
        jIndex
        tensor
        microName
        tensorCase
        tensName
    end

    methods (Access = public)

        function obj = VademecumLoader(d)
            obj.init(d)
            obj.loadMxMy();
            obj.loadHomogenizedTensor()
            obj.loadAmplificatorTensor()
        end

    end

    methods (Access = private)

        function init(obj,d)
            obj.index = d.index;
            obj.path = d.path;
            obj.microName = d.microName;
        end

        function loadMxMy(obj)
            m1Path = fullfile(obj.path,'m1m2.txt');
            obj.dataBase.mx = load(m1Path);
            obj.dataBase.my = obj.dataBase.mx;
            obj.dataBase.nmx = length(obj.dataBase.mx);
            obj.dataBase.nmy = length(obj.dataBase.my);
        end

        function loadHomogenizedTensor(obj)
            obj.tensorCase = 'C_';
            obj.tensName = 'ConstitutiveTensor';
            obj.loadTensors();
            obj.dataBase.C = obj.tensor;
        end

        function loadAmplificatorTensor(obj)
            obj.tensorCase = 'P_';
            obj.tensName = 'AmplificatorTensor';
            obj.loadTensors();
            obj.dataBase.P = obj.tensor;
        end

        function loadTensors(obj)
            nPlot = length(obj.index);
            for iplot = 1:nPlot
                obj.obtainIJindex(iplot);
                obj.loadTensor();
            end
        end

       function obtainIJindex(obj,iplot)
            obj.iIndex = obj.index(iplot,1);
            obj.jIndex = obj.index(iplot,2);
       end

        function loadTensor(obj)
            mN    = obj.microName;
            tName = obj.tensName;
            i = obj.iIndex;
            j = obj.jIndex;
            name = [obj.tensorCase,num2str(i),num2str(j)];
            fullPath = fullfile(obj.path,mN,tName,[mN,name,'.txt']);
            obj.tensor(i,j,:,:) = load(fullPath);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan