text:
using Test
using Base.BinaryPlatforms
using BinaryBuilderBase
using BinaryBuilderBase: abi_agnostic, get_concrete_platform, march, platform_dlext, platform_exeext,
                         nbits, proc_family

@testset "Supported Platforms" begin
    all = supported_platforms()
    opt_out_specific = supported_platforms(exclude=[Platform("x86_64", "linux"; libc="glibc")])
    opt_out_fx = supported_platforms(exclude=Sys.islinux)

    @test length(all) == length(opt_out_specific)+1
    @test any(opt_out_specific .== [Platform("i686", "linux", libc="glibc")])
    @test !any(opt_out_fx .== [Platform("i686", "linux", libc="glibc")])

    @test sort(
        [
            Platform("x86_64", "windows"),
            Platform("i686", "linux"; libc="musl"),
            Platform("i686", "linux"; libc="glibc"),
            Platform("x86_64", "macos")
        ], by = triplet) == [
            Platform("i686", "linux"; libc="glibc"),
            Platform("i686", "linux"; libc="musl"),
            Platform("x86_64", "macos"),
            Platform("x86_64", "windows"),
        ]
end

@testset "abi_agnostic" begin
    p = Platform("x86_64", "linux"; libc="musl", libgfortran_version=v"3", cxxstring_abi="cxx03", libstdcxx_version=v"3.4.26")
    @test abi_agnostic(p) == Platform("x86_64", "linux"; libc="musl")
end

@testset "AnyPlatform" begin
    # Test some AnyPlatform properties
    @test triplet(AnyPlatform()) == "any"
    @test abi_agnostic(AnyPlatform()) === AnyPlatform()
    @test arch(AnyPlatform()) == "any"
    @test repr(AnyPlatform()) == "AnyPlatform"

    # In the build environment we want AnyPlatform to look like x86_64-linux-musl
    @test get_concrete_platform(
            AnyPlatform();
            compilers = [:c],
            preferred_gcc_version = v"7",
            preferred_llvm_version = v"9",
        ) == get_concrete_platform(
            Platform("x86_64", "linux"; libc="musl");
            compilers = [:c],
            preferred_gcc_version = v"7",
            preferred_llvm_version = v"9",
        )
    @test BinaryBuilderBase.choose_shards(AnyPlatform()) == BinaryBuilderBase.choose_shards(Platform("x86_64", "linux"; libc="musl"))
    @test BinaryBuilderBase.aatriplet(AnyPlatform()) == BinaryBuilderBase.aatriplet(Platform("x86_64", "linux"; libc="musl"))
end

@testset "Target properties" begin
    for p in [Platform("i686", "linux"), Platform("i686", "windows"), Platform("armv7l", "linux")]
        @test nbits(p) == 32
    end

    for p in [Platform("x86_64", "linux"), Platform("x86_64", "windows"), Platform("aarch64", "linux"),
              Platform("powerpc64le", "linux"), Platform("x86_64", "macos")]
        @test nbits(p) == 64
    end

    for p in [Platform("x86_64", "linux"), Platform("x86_64", "macos"), Platform("i686", "windows")]
        @test proc_family(p) == "intel"
    end
    for p in [Platform("aarch64", "linux"; libc="musl"), Platform("armv7l", "linux")]
        @test proc_family(p) == "arm"
    end
    @test proc_family(Platform("powerpc64le", "linux")) == "power"

    for p in [Platform("x86_64", "linux"), Platform("x86_64", "freebsd"), Platform("powerpc64le", "linux"), Platform("x86_64", "macos")]
        @test platform_exeext(p) == ""
    end
    @test platform_exeext(Platform("x86_64", "windows")) == ".exe"
    @test platform_exeext(Platform("i686", "windows")) == ".exe"

    @test platform_dlext(Platform("aarch64", "linux"; march="thunderx2")) == "so"
    @test platform_dlext(Platform("x86_64", "freebsd"; march="x86_64")) == "so"
    @test platform_dlext(Platform("x86_64", "macos"; march="avx512")) == "dylib"
    @test platform_dlext(Platform("i686", "windows"; cuda="10.1")) == "dll"
end

using BinaryBuilderBase: get_march_flags, get_all_arch_names, get_all_march_names
@testset "march flags" begin
    # test one that is common between gcc and clang
    @test get_march_flags("x86_64", "avx", "gcc") == ["-march=sandybridge", "-mtune=sandybridge"]

    # Test one that is different between gcc and clang
    @test get_march_flags("aarch64", "armv8_2_crypto", "gcc") == ["-march=armv8-a+lse+crc+rdma+aes+sha2", "-mtune=cortex-a76"]
    @test get_march_flags("aarch64", "armv8_2_crypto", "clang") == ["-march=armv8-a+lse+crc+rdm+aes+sha2", "-mtune=cortex-a76"]

    @test get_march_flags("aarch64", nothing, "gcc") == get_march_flags("aarch64", "armv8_0", "gcc")
    @test get_march_flags("armv7l", nothing, "clang") == get_march_flags("armv7l", "armv7l", "clang")
    @test get_march_flags("x86_64", nothing, "clang") == get_march_flags("x86_64", "x86_64", "clang")

    # Get all architectures and all microarchitectures for the different architectures
    @test sort(get_all_arch_names()) == ["aarch64", "armv6l", "armv7l", "i686", "powerpc64le", "x86_64"]
    @test sort(get_all_march_names("x86_64")) == ["avx", "avx2", "avx512", "x86_64"]
    @test sort(get_all_march_names("armv7l")) == ["armv7l", "neonvfpv4"]
end


text:
function CubedGrid(n,OrientFace,Param)
  Grid = GridStruct(;)
  Grid.nBar=[ 0  1   0   1
             -1  0  -1   0];
  Grid.Dim=3;
  Grid.Type="Quad";
  Grid.Rad=Param.RadEarth;
  Grid.Form="Sphere";
  dd=2.0e0/n;
  NumNodes=(6*(n-1)*(n-1)+12*(n-1)+8);
  Nodes = map(1:NumNodes) do i
    Node([0 0 0],0);
  end
  NodeNumber=1;
  #Faces
  # West
  NodeNumberW=NodeNumber;
  x=-1.0;
  y=0.0;
  z=0.0;
  for k=1:n-1
    for j=1:n-1
      Nodes[NodeNumber]=Node(CubePoint(-1,j,k,n,x,y,z,Grid.Rad),NodeNumber);
      NodeNumber=NodeNumber+1;
    end
  end
  # ! East
  NodeNumberE=NodeNumber;
  x=1.0;
  y=0.0;
  z=0.0;
  for k=1:n-1
    for j=1:n-1
      Nodes[NodeNumber]=Node(CubePoint(-1,j,k,n,x,y,z,Grid.Rad),NodeNumber);
      NodeNumber=NodeNumber+1;
    end
  end
  # ! South
  NodeNumberS=NodeNumber;
  x=0.0;
  y=-1.0;
  z=0.0;
  for k=1:n-1
    for i=1:n-1
      Nodes[NodeNumber]=Node(CubePoint(i,-1,k,n,x,y,z,Grid.Rad),NodeNumber);
      NodeNumber=NodeNumber+1;
    end
  end
  # ! North
  NodeNumberN=NodeNumber;
  x=0.0;
  y=1.0;
  z=0.0;
  for k=1:n-1
    for i=1:n-1
      Nodes[NodeNumber]=Node(CubePoint(i,-1,k,n,x,y,z,Grid.Rad),NodeNumber);
      NodeNumber=NodeNumber+1;
    end
  end
  # ! Bottom
  NodeNumberB=NodeNumber;
  x=0.0;
  y=0.0;
  z=-1.0;
  for j=1:n-1
    for i=1:n-1
      Nodes[NodeNumber]=Node(CubePoint(i,j,-1,n,x,y,z,Grid.Rad),NodeNumber);
      NodeNumber=NodeNumber+1;
    end
  end
  # ! Top
  NodeNumberT=NodeNumber;
  x=0.0;
  y=0.0;
  z=1.0;
  for j=1:n-1
    for i=1:n-1
      Nodes[NodeNumber]=Node(CubePoint(i,j,-1,n,x,y,z,Grid.Rad),NodeNumber);
      NodeNumber=NodeNumber+1;
    end
  end
  # !Edges
  # !West East
  NodeNumberWEmm=NodeNumber;
  x=0.0;
  y=-1.0;
  z=-1.0;
  for i=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(i,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end
  NodeNumberWEpm=NodeNumber;
  x=0.0;
  y=1.0;
  z=-1.0;
  for i=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(i,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end
  NodeNumberWEmp=NodeNumber;
  x=0.0;
  y=-1.0;
  z=1.0;
  for i=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(i,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end
  NodeNumberWEpp=NodeNumber;
  x=0.0;
  y=1.0;
  z=1.0;
  for i=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(i,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end
  # South North
  NodeNumberSNmm=NodeNumber;
  x=-1.0;
  y=0.0;
  z=-1.0;
  for j=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(-1,j,-1,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end
  NodeNumberSNpm=NodeNumber;
  x=1.0;
  y=0.0;
  z=-1.0;
  for j=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(-1,j,-1,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end
  NodeNumberSNmp=NodeNumber;
  x=-1.0;
  y=0.0;
  z=1.0;
  for j=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(-1,j,-1,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end
  NodeNumberSNpp=NodeNumber;
  x=1.0;
  y=0.0;
  z=1.0;
  for j=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(-1,j,-1,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end
  # Bottom Top
   NodeNumberBTmm=NodeNumber;
   x=-1.0;
  y=-1.0;
  z=0.0;
  for k=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(-1,-1,k,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end

   NodeNumberBTpm=NodeNumber;
    x=1.0;
  y=-1.0;
  z=0.0;
  for k=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(-1,-1,k,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end

   NodeNumberBTmp=NodeNumber;
    x=-1.0;
  y=1.0;
  z=0.0;
  for k=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(-1,-1,k,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end

  NodeNumberBTpp=NodeNumber;
   x=1.0;
  y=1.0;
  z=0.0;
  for k=1:n-1
    Nodes[NodeNumber]=Node(CubePoint(-1,-1,k,n,x,y,z,Grid.Rad),NodeNumber);
    NodeNumber=NodeNumber+1;
  end
  # Nodes
  NodeNumbermmm=NodeNumber;
  x=-1.0e0;
  y=-1.0e0;
  z=-1.0e0;
  Nodes[NodeNumber]=Node(CubePoint(-1,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
  NodeNumber=NodeNumber+1;
  NodeNumberpmm=NodeNumber;
  x=1.0e0;
  y=-1.0e0;
  z=-1.0e0;
  Nodes[NodeNumber]=Node(CubePoint(-1,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
  NodeNumber=NodeNumber+1;
  NodeNumbermpm=NodeNumber;
  x=-1.0e0;
  y=1.0e0;
  z=-1.0e0;
  Nodes[NodeNumber]=Node(CubePoint(-1,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
  NodeNumber=NodeNumber+1;
  NodeNumberppm=NodeNumber;
  x=1.0e0;
  y=1.0e0;
  z=-1.0e0;
  Nodes[NodeNumber]=Node(CubePoint(-1,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
  NodeNumber=NodeNumber+1;
  NodeNumbermmp=NodeNumber;
  x=-1.0e0;
  y=-1.0e0;
  z=1.0e0;
  Nodes[NodeNumber]=Node(CubePoint(-1,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
  NodeNumber=NodeNumber+1;
  NodeNumberpmp=NodeNumber;
  x=1.0e0;
  y=-1.0e0;
  z=1.0e0;
  Nodes[NodeNumber]=Node(CubePoint(-1,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
  NodeNumber=NodeNumber+1;

  NodeNumbermpp=NodeNumber;
  x=-1.0e0;
  y=1.0e0;
  z=1.0e0;
  Nodes[NodeNumber]=Node(CubePoint(-1,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
  NodeNumber=NodeNumber+1;
  NodeNumberppp=NodeNumber;
  x=1.0e0;
  y=1.0e0;
  z=1.0e0;
  Nodes[NodeNumber]=Node(CubePoint(-1,-1,-1,n,x,y,z,Grid.Rad),NodeNumber);
  NodeNumber=NodeNumber+1;
  Grid.Nodes=Nodes;

  NumEdges=12*(n-1)*n+12*n;
  NumEdgesI=12*(n-1)*n+12*n;
  Edges = map(1:NumEdges) do i
    Edge([1 2],Grid,0,0,"",0);
  end
  EdgeNumber=1;
  # West
  (Edges,EdgeNumber,EdgeNumberW1,EdgeNumberW2)=InsertFaceEdge(n,EdgeNumber,NodeNumberW
    ,NodeNumberBTmm,NodeNumberBTmp
    ,NodeNumberSNmm,NodeNumberSNmp
    ,Edges,Grid);

  # East
  (Edges,EdgeNumber,EdgeNumberE1,EdgeNumberE2)=InsertFaceEdge(n,EdgeNumber,NodeNumberE
                     ,NodeNumberBTpm,NodeNumberBTpp
                     ,NodeNumberSNpm,NodeNumberSNpp
                     ,Edges,Grid);
  # South
  (Edges,EdgeNumber,EdgeNumberS1,EdgeNumberS2)=InsertFaceEdge(n,EdgeNumber,NodeNumberS
                     ,NodeNumberBTmm,NodeNumberBTpm
                     ,NodeNumberWEmm,NodeNumberWEmp
                     ,Edges,Grid);
  # North
  (Edges,EdgeNumber,EdgeNumberN1,EdgeNumberN2)=InsertFaceEdge(n,EdgeNumber,NodeNumberN
                     ,NodeNumberBTmp,NodeNumberBTpp
                     ,NodeNumberWEpm,NodeNumberWEpp
                     ,Edges,Grid);
  # Bottom
  (Edges,EdgeNumber,EdgeNumberB1,EdgeNumberB2)=InsertFaceEdge(n,EdgeNumber,NodeNumberB
                     ,NodeNumberSNmm,NodeNumberSNpm
                     ,NodeNumberWEmm,NodeNumberWEpm
                     ,Edges,Grid);
   # Top
   (Edges,EdgeNumber,EdgeNumberT1,EdgeNumberT2)=InsertFaceEdge(n,EdgeNumber,NodeNumberT
     ,NodeNumberSNmp,NodeNumberSNpp
     ,NodeNumberWEmp,NodeNumberWEpp
     ,Edges,Grid);

# Edges
#West East
  (Edges,EdgeNumber,EdgeNumberWEmm)=InsertEdgeEdge(n,EdgeNumber,NodeNumberWEmm,NodeNumbermmm,NodeNumberpmm,Edges,Grid);
  (Edges,EdgeNumber,EdgeNumberWEpm)=InsertEdgeEdge(n,EdgeNumber,NodeNumberWEpm,NodeNumbermpm,NodeNumberppm,Edges,Grid);
  (Edges,EdgeNumber,EdgeNumberWEmp)=InsertEdgeEdge(n,EdgeNumber,NodeNumberWEmp,NodeNumbermmp,NodeNumberpmp,Edges,Grid);
  (Edges,EdgeNumber,EdgeNumberWEpp)=InsertEdgeEdge(n,EdgeNumber,NodeNumberWEpp,NodeNumbermpp,NodeNumberppp,Edges,Grid);
# South North
  (Edges,EdgeNumber,EdgeNumberSNmm)=InsertEdgeEdge(n,EdgeNumber,NodeNumberSNmm,NodeNumbermmm,NodeNumbermpm,Edges,Grid);
  (Edges,EdgeNumber,EdgeNumberSNpm)=InsertEdgeEdge(n,EdgeNumber,NodeNumberSNpm,NodeNumberpmm,NodeNumberppm,Edges,Grid);
  (Edges,EdgeNumber,EdgeNumberSNmp)=InsertEdgeEdge(n,EdgeNumber,NodeNumberSNmp,NodeNumbermmp,NodeNumbermpp,Edges,Grid);
  (Edges,EdgeNumber,EdgeNumberSNpp)=InsertEdgeEdge(n,EdgeNumber,NodeNumberSNpp,NodeNumberpmp,NodeNumberppp,Edges,Grid);
# Bottom Top
  (Edges,EdgeNumber,EdgeNumberBTmm)=InsertEdgeEdge(n,EdgeNumber,NodeNumberBTmm,NodeNumbermmm,NodeNumbermmp,Edges,Grid);
  (Edges,EdgeNumber,EdgeNumberBTpm)=InsertEdgeEdge(n,EdgeNumber,NodeNumberBTpm,NodeNumberpmm,NodeNumberpmp,Edges,Grid);
  (Edges,EdgeNumber,EdgeNumberBTmp)=InsertEdgeEdge(n,EdgeNumber,NodeNumberBTmp,NodeNumbermpm,NodeNumbermpp,Edges,Grid);
  (Edges,EdgeNumber,EdgeNumberBTpp)=InsertEdgeEdge(n,EdgeNumber,NodeNumberBTpp,NodeNumberppm,NodeNumberppp,Edges,Grid);

  Grid.Edges=Edges;

  NumFaces=6*n*n;
  Faces = map(1:NumFaces) do i
    first(Face([0 0 0 0],Grid,0,"x"));
  end

  FaceNumber=1;
# Faces
# West
  (Faces,FaceNumber,Grid)=InsertFaceFace(n,FaceNumber,"W",EdgeNumberW1,EdgeNumberW2
                     ,EdgeNumberSNmm,EdgeNumberSNmp
                     ,EdgeNumberBTmm,EdgeNumberBTmp,Grid,Faces,OrientFace);

# East
  (Faces,FaceNumber,Grid)=InsertFaceFace(n,FaceNumber,"E",EdgeNumberE1,EdgeNumberE2
                     ,EdgeNumberSNpm,EdgeNumberSNpp
                     ,EdgeNumberBTpm,EdgeNumberBTpp,Grid,Faces,OrientFace);
# South
  (Faces,FaceNumber,Grid)=InsertFaceFace(n,FaceNumber,"S",EdgeNumberS1,EdgeNumberS2
                     ,EdgeNumberWEmm,EdgeNumberWEmp
                     ,EdgeNumberBTmm,EdgeNumberBTpm,Grid,Faces,OrientFace);
# North
  (Faces,FaceNumber,Grid)=InsertFaceFace(n,FaceNumber,"N",EdgeNumberN1,EdgeNumberN2
                     ,EdgeNumberWEpm,EdgeNumberWEpp
                     ,EdgeNumberBTmp,EdgeNumberBTpp,Grid,Faces,OrientFace);
# Bottom
  (Faces,FaceNumber,Grid)=InsertFaceFace(n,FaceNumber,"B",EdgeNumberB1,EdgeNumberB2
                     ,EdgeNumberWEmm,EdgeNumberWEpm
                     ,EdgeNumberSNmm,EdgeNumberSNpm,Grid,Faces,OrientFace);
# Top
  (Faces,FaceNumber,Grid)=InsertFaceFace(n,FaceNumber,"T",EdgeNumberT1,EdgeNumberT2
                     ,EdgeNumberWEmp,EdgeNumberWEpp
                     ,EdgeNumberSNmp,EdgeNumberSNpp,Grid,Faces,OrientFace);

  Grid.Faces=Faces;

  Grid.NumNodes=size(Grid.Nodes,1);
  Grid.NumEdges=size(Grid.Edges,1);
  Grid.NumEdgesI=size(Grid.Edges,1);
  Grid.NumEdgesB=0;
  Grid.NumFaces=size(Grid.Faces,1);
  Grid.Dim=3;
  Grid=Orientation(Grid);
  Grid=Renumbering(Grid);
  Grid=FacesInNodes(Grid);

end

function CubePoint(i1,i2,i3,n,x,y,z,Rad)

if i1>0
  x=tan(i1*pi/(2*n)-0.25*pi);
end
if i2>0
  y=tan(i2*pi/(2*n)-0.25*pi);
end
if i3>0
  z=tan(i3*pi/(2*n)-0.25*pi);
end
N=[x;y;z];
N=N/norm(N)*Rad;
return N
end

function InsertFaceEdge(n,EdgeNumber,NodeNumberStart,
                         NodeNumberE1Start1,NodeNumberE2Start1,
                         NodeNumberE1Start2,NodeNumberE2Start2,Edges,Grid)

  NodeNumber=NodeNumberStart;
  NodeNumberE1=NodeNumberE1Start1;
  NodeNumberE2=NodeNumberE2Start1;
  EdgeNumberStart1=EdgeNumber;
  for j=1:n-1
    for i=1:n
      N1=NodeNumber-1;
      N2=NodeNumber;
      if i==1
        N1=NodeNumberE1;
        NodeNumberE1=NodeNumberE1+1;
      end
      if i==n
        N2=NodeNumberE2;
        NodeNumberE2=NodeNumberE2+1;
        NodeNumber=NodeNumber-1;
      end
      Edges[EdgeNumber]=Edge([N1 N2],Grid,EdgeNumber,EdgeNumber,"X",EdgeNumber);
      EdgeNumber=EdgeNumber+1;
      NodeNumber=NodeNumber+1;
    end
  end
  NodeNumber=NodeNumberStart;
  NodeNumberE1=NodeNumberE1Start2;
  NodeNumberE2=NodeNumberE2Start2;
  EdgeNumberStart2=EdgeNumber;
  for j=1:n
    for i=1:n-1
      N1=NodeNumber-(n-1);
      N2=NodeNumber;
      if j==1
        N1=NodeNumberE1;
        NodeNumberE1=NodeNumberE1+1;
      end
      if j==n
        N2=NodeNumberE2;
        NodeNumberE2=NodeNumberE2+1;
      end
      Edges[EdgeNumber]=Edge([N1 N2],Grid,EdgeNumber,EdgeNumber,"X",EdgeNumber);
      EdgeNumber=EdgeNumber+1;
      NodeNumber=NodeNumber+1;
    end
  end
  return (Edges,EdgeNumber,EdgeNumberStart1,EdgeNumberStart2)
end

function InsertEdgeEdge(n,EdgeNumber,
    NodeNumberStart,NodeNumberE1,NodeNumberE2,
     Edges,Grid)
  NodeNumber=NodeNumberStart;
  EdgeNumberStart=EdgeNumber;
  for i=1:n
    N1=NodeNumber-1;
    N2=NodeNumber;
    if i==1
      N1=NodeNumberE1;
    end
    if i==n
      N2=NodeNumberE2;
    end
    Edges[EdgeNumber]=Edge([N1 N2],Grid,EdgeNumber,EdgeNumber,"X",EdgeNumber);
    EdgeNumber=EdgeNumber+1;
    NodeNumber=NodeNumber+1;
  end
  return (Edges,EdgeNumber,EdgeNumberStart)
end

function InsertFaceFace(n,FaceNumber,Type,EdgeNumberStart1,EdgeNumberStart2,
                         EdgeNumberStartEW1,EdgeNumberStartEW2,
                         EdgeNumberStartSN1,EdgeNumberStartSN2,
                         Grid,Faces,OrientFace)

  EdgeNumber1=EdgeNumberStart1;
  EdgeNumber2=EdgeNumberStart2;
  EdgeNumberEW1=EdgeNumberStartEW1;
  EdgeNumberEW2=EdgeNumberStartEW2;
  EdgeNumberSN1=EdgeNumberStartSN1;
  EdgeNumberSN2=EdgeNumberStartSN2;
  for j=1:n
    for i=1:n
      if (i==1 && j==1) || (i==1 && j==n)||
          (i==n && j==1) || (i==n && j==n)

      end
      #Face(FaceNumber)#Number=FaceNumber
      #Face(FaceNumber)#Boundary=1
      #ALLOCATE(Face(FaceNumber)#Edge(4))
      E1=EdgeNumber1-n;
      E2=EdgeNumber1;
      E3=EdgeNumber2-1;
      E4=EdgeNumber2;
      if j==1
        E1=EdgeNumberEW1;
        EdgeNumberEW1=EdgeNumberEW1+1;
      end
      if j==n
        E2=EdgeNumberEW2;
        EdgeNumberEW2=EdgeNumberEW2+1;
      end
      if i==1
        E3=EdgeNumberSN1;
        EdgeNumberSN1=EdgeNumberSN1+1;
      end
      if i==n
        E4=EdgeNumberSN2;
        EdgeNumberSN2=EdgeNumberSN2+1;
        EdgeNumber2=EdgeNumber2-1;
      end
      EdgeNumber1=EdgeNumber1+1;
      EdgeNumber2=EdgeNumber2+1;
      (Faces[FaceNumber],Grid)=Face([E1 E2 E3 E4],Grid,FaceNumber,Type,OrientFace);
      FaceNumber=FaceNumber+1;
    end
  end
  return (Faces,FaceNumber,Grid)
end

ext: jl


text:
cdir = pwd(); ii = findlast("MDP", cdir); MDPpath = cdir[1:ii[end]] * "/";
include(MDPpath * "Installation/setup.jl");

# Read xyz file to obtain configurations
filename = "/Users/ngoccuongnguyen/Dropbox (MIT)/MDP/data/LJCluster/curated_lj_cluster.xyz";
species = ["Ar"];
config = readEXTXYZ(filename, species);

# shift atom positions to put them in the box [0,10]x[0,10]x[0,10]
config.x = config.x .+ 5.0

# cumalative sum of numbers of atoms
natom = [0; cumsum(config.natom[:])];

ci = 1; # choose a configuration to start MD simulation
x = config.x[:,(natom[ci]+1):natom[ci+1]] # atom positions of configuration ci
t = config.t[:,(natom[ci]+1):natom[ci+1]] # atom types of configuration ci
q = [];                                   # atom charges of configuration ci
v = 0.0*x # zero velocity

# 3 vectors define the simulation box
a = [10.0; 0.0; 0.0]
b = [0.0; 10.0; 0.0]
c = [0.0; 0.0; 10.0]

# periodic boundary conditions
pbc = [0; 0; 0]

# Bessel parameters
mu = [9080.776024081284, 22740.90481568728, 22826.411376542143, 12272.607095393349, 3542.9517448941156, 437.16635703260755]
rcut = 3.294988099200801
alpha = 4.902241343306524

# define Bessel potential
potentials[1] = addpotential("pair", "nonbonded", [0], "scaledbessel", mu, rcut);
include(potentials[1].potentialfunction * ".jl");  # include the potentials file

# SNAP parameters
twojmax = 4
ntypes = 1
rcutfac = rcut;
rfac0 = 1.0;
rmin0 = 0;
bzeroflag = 0;
quadraticflag = 0;
switchflag = 1;
chemflag = 0;
bnormflag = 0;
wselfallflag = 0;

snaparam = [ntypes, twojmax, rcutfac, rfac0, rmin0, bzeroflag, switchflag, quadraticflag, chemflag, bnormflag, wselfallflag]
elemradius = 0.5*ones(Float64, ntypes)
elemweight = ones(Float64, ntypes)

snapcoeff = [0.0
9.51697531984316e-5
 -3.503209826706362e-6
 -0.00015477632687816527
  8.864246681893688e-5
 -6.228723743014278e-5
  0.0002572563742710612
  1.0535143289316965e-6
  9.529421498586589e-5
  9.123169496032785e-6
 -0.0001693491133185178
  0.00013696751121984198
  3.328433092023331e-5
 -6.500144488870558e-5
 -3.968313413722244e-5]

# define SNAP potential
potentials[2] = initsna(snaparam, elemradius, elemweight, snapcoeff)

app.unitstyle = "lj"; # unit system
app.ensemblemode = "nve";
app.time = 0.0;             # initial time
app.dt = 0.5e-3;            # time step size
app.ntimesteps = 1000;       # time steps
app.globalfreq = 10;        # time frequency to print out global outputs
app.neighskin = 0.5;        # add a skin distance for neighborlist
app.eta = [rcut alpha]
app.atommasses = [1.0 1.0];

# run MD simulation using velocity Verlet algorithm
VelocityVerlet(x, v, q, t, a, b, c, pbc, app, potentials)

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
export jlboost!, jlboost

using DataFrames: nrow, ncol
using Tables

using ..JLBoostTrees: JLBoostTree

"""
    jlboost(df, target, features = setdiff(names(df), (target, prev_w, new_weight)),
        warm_start = fill(0.0, nrow(df)); nrounds = 1, eta = 0.3, lambda = 0, gamma = 0,
        max_depth = 6, subsample = 1, colsample_bytree=1, colsample_bylevel=1, colsample_bynode=1,
        verbose = false)

Fit a tree boosting model with a DataFrame, df, and target symbol and allowed features.

This is based on the xgboost interface, where possible the parameters have the same name as xgboost,
see https://xgboost.readthedocs.io/en/latest/parameter.html

* nrounds: Number of trees to fit
* warmstart: A vector of weights from which to start training. Defaults to 0. The warmstart may be
    different for every row. This is designed to allow the model to improve upon existing models.
* eta: The learning rate. Also known as the weight of each tree in the final summation of trees
* lambda: XGBoost lambda hyperparameter
* gamma: XGBoost gamma hyperparameter
* max_depth: the maximum depth of each tree
* subsample: 0-1, the proportion of rows to subsample for each tree build
* verbose: Print more information
* colsample_bytree: (0-1] The proportion of feature column to sample for each tree.
* min_child_weight: The weight that needs to be in each child node before a split can occur. The
    weight is the hessian (2nd derivative) of the loss function, which happens to be 1 for squares
    loss.
* colsample_bylevel: Not yet implemented
* colsample_bynode: Not yet implemented
* monotone_contraints: Not yet implemented
* interaction_constraints: Not yet implemented
"""
function jlboost(df, target::Union{Symbol, String}; kwargs...)
    target = Symbol(target)
    warm_start = fill(0.0, nrow(df))
	jlboost(df, target, setdiff(Tables.columnnames(df), [target]), warm_start; kwargs...)
end

function jlboost(df, target::Union{Symbol, String}, warm_start::AbstractVector{T}; kwargs...) where T <: Number
    target = Symbol(target)
	jlboost(df, target, setdiff(names(df), [target]), warm_start)
end

function jlboost(df, target::Union{Symbol, String}, features::AbstractVector{T}; kwargs...) where T <: Union{String, Symbol}
    target = Symbol(target)
    features = Symbol.(features)
	jlboost(df, target, features, fill(0.0, nrow(df)); kwargs...)
end

function jlboost(df, target::Union{Symbol, String}, features::AbstractVector,
    warm_start::AbstractVector, loss = LogitLogLoss();
    subsample = 1, colsample_bytree = 1, max_depth = 6, max_leaves = 0, kwargs...)

    @assert 0 < subsample <= 1
    @assert 0 < colsample_bytree <= 1
    @assert Tables.istable(df)

    target = Symbol(target)
    features = Symbol.(features)

    # a sample of the rows
    row_sampling_bytree_strategy = select_row_sampling_strategy(subsample)

    # a function to sample the columns
    col_sampling_bytree_strategy = select_col_sampling_strategy(colsample_bytree)

    if max_leaves > 0
        if max_depth > 0
            @warn "You have set max_leaves=$max_leaves but max_depth > 0. The max_depth parameter is ignored."
        end
        tree_growth = lossguide
        stopping_criterion = max_leaves_stopping_criterion(max_leaves)
    else
        tree_growth = depth_wise
        stopping_criterion = max_depth_stopping_criterion(max_depth)
    end

    # TODO look at target column and provide a possible selection of loss
    # e.g. if the target is numeric then RSMELoss is more appropriate
    jlboost(df, target, features, warm_start, loss,
            row_sampling_bytree_strategy,
            col_sampling_bytree_strategy,
            tree_growth,
            stopping_criterion; kwargs...)
end

# the most canonical version of jlboost is here
function jlboost(df, target, features, warm_start::AbstractVector,
    loss,
    row_sampling_strategy::Function,
    col_sampling_bytree_strategy::Function,
    tree_growth::Function,
    stopping_criterion::Function;
	nrounds = 1, eta = 1.0, verbose = false, kwargs...)
    # eta = 1, lambda = 0, gamma = 0,  min_child_weight = 1, colsample_bylevel = 1, colsample_bynode = 1,
	#, ,  colsample_bynode = 1,

    @assert nrounds >= 1
	@assert Tables.istable(df)

    target = Symbol(target)
    features = Symbol.(features)

	dfc = Tables.columns(df)

    # res_jlt = result JLBoost trees
	res_jlt = AbstractJLBoostTree[]

    # fit the next round
	for nround in 1:nrounds
		if verbose
			println("Fitting tree #$(nround)")
		end

        # sample new columns
		features_sample = col_sampling_bytree_strategy(features, df, target, warm_start, loss;
                                                   nrounds=nrounds, eta=eta, kwargs...)
        # dfs = DataFrame Sampled
        dfs = row_sampling_strategy(dfc)
        if nround == 1
            warm_start = fill(0.0, nrow(dfs))
        else
            warm_start = predict(res_jlt[1:nround-1], dfs)
        end

        println(nround)
        println(dfc)

        new_jlt = _fit_tree!(loss, dfc, target, features_sample, warm_start, JLBoostTree(0.0),
                             tree_growth,
                             stopping_criterion; verbose=verbose, kwargs...);

        println("mehmehmehmeh")
        # added a new round of tree
        push!(res_jlt, eta*deepcopy(new_jlt))
	end
	res_jlt

    JLBoostTreeModel(res_jlt, loss, target)
end


text:
#------------------------------------------------------------------------------
"""
    printmat([fh::IO],x...;colNames=[],rowNames=[],
             width=10,prec=3,NoPrinting=false,StringFmt="",cell00="")

Print all elements of a matrix (or several) with predefined formatting. It can also handle
OffsetArrays. StringFmt = "csv" prints using a csv format.

# Input
- `fh::IO`:            (optional) file handle. If not supplied, prints to screen
- `x::Array(s)`:       (of numbers, dates, strings, ...) to print
- `colNames::Array`:   of strings with column headers
- `rowNames::Array`:   of strings with row labels
- `width::Int`:        (keyword) scalar, minimum width of printed cells
- `prec::Int`:         (keyword) scalar, precision of printed cells
- `NoPrinting::Bool`:  (keyword) bool, true: no printing, just return formatted string [false]
- `StringFmt::String`: (keyword) string, "", "csv"
- `cell00::String`:    (keyword) string, for row 0, column 0

# Output
- str         (if NoPrinting) string, (otherwise nothing)

# Examples
```
x = [11 12;21 22]
printmat(x)
```
```
x = [1 "ab"; Date(2018,10,7) 3.14]
printmat(x,width=20,colNames=["col 1","col 2"])
```
```
printmat([11,12],[21,22])
```
Can also call as
```
opt = Dict(:rowNames=>["1";"4"],:width=>10,:prec=>3,:NoPrinting=>false,:StringFmt=>"")
printmat(x;colNames=["a","b"],opt...)     #notice ; and ...
```
(not all keywords are needed)

# Requires
- fmtNumPs

# Notice
- The prefixN and suffixN could potentially be made function inputs. This would allow
a fairly flexible way to format tables.

Paul.Soderlind@unisg.ch

"""
function printmat(fh::IO,x...;colNames=[],rowNames=[],
                  width=10,prec=3,NoPrinting=false,StringFmt="",cell00="")

  isempty(x) && return nothing                         #do nothing is isempty(x)

  typeTestQ = any(!=(eltype(x[1])),[eltype(z) for z in x])  #test if eltype(x[i]) differs
  if typeTestQ                                      #create matrix from tuple created by x...
    x = hcat(Matrix{Any}(hcat(x[1])),x[2:end]...)   #preserving types of x[i]
  else
    x = hcat(x...)
  end

  (m,n) = (size(x,1),size(x,2))

  (length(rowNames) == 1 < m) && (rowNames = [string(rowNames[1],i) for i = 1:m])  #"ri"
  (length(colNames) == 1 < n) && (colNames = [string(colNames[1],i) for i = 1:n])  #"ci"

  if StringFmt == "csv"
    (prefixN,suffixN)   = (fill("",n),vcat(fill(",",n-1),""))  #prefix and suffix for column 1:n
    (prefixC0,suffixC0) = ("",",")                             #prefix and suffix for column 0
  else
    (prefixN,suffixN) = (fill("",n),fill("",n))
    (prefixC0,suffixC0) = ("","")
  end

  if length(rowNames) == 0                         #width of column 0 (cell00 and rowNames)
    col0Width = 0
  else
    col0Width = maximum(length,vcat(cell00,rowNames)) + length(prefixC0) + length(suffixC0)
  end

  colWidth = [width + length(prefixN[j]) + length(suffixN[j]) for j=1:n]  #widths of column 1:n

  iob = IOBuffer()

  if !isempty(colNames)                                #print (cell00,colNames), if any
    !isempty(cell00) ?  txt0 = string(prefixC0,cell00,suffixC0) : txt0 = ""
    print(iob,rpad(txt0,col0Width))
    for j = 1:n                                #loop over columns
      print(iob,lpad(string(prefixN[j],colNames[j],suffixN[j]),colWidth[j]))
    end
    print(iob,"\n")
  end
                                                       #print rowNames and x
  (i0,j0) = (1 - first(axes(x,1)),1 - first(axes(x,2)))   #i+i0,j+j0 give traditional indices
  for i in axes(x,1)                           #loop over rows
    !isempty(rowNames) && print(iob,rpad(string(prefixC0,rowNames[i+i0],suffixC0),col0Width))
    for j in axes(x,2)                         #loop over columns
      print(iob,fmtNumPs(x[i,j],width,prec,"right",prefix=prefixN[j+j0],suffix=suffixN[j+j0]))
    end
    print(iob,"\n")
  end
  str = String(take!(iob))

  if NoPrinting                              #no printing, just return str
    return str
  else                                       #print, return nothing
    print(fh,str,"\n")
    return nothing
  end

end
                        #when fh is not supplied: printing to screen
printmat(x...;colNames=[],rowNames=[],width=10,prec=3,NoPrinting=false,StringFmt="",cell00="") =
    printmat(stdout::IO,x...;colNames,rowNames,width,prec,NoPrinting,StringFmt,cell00)
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    printlnPs([fh::IO],z...;width=10,prec=3)

Subsitute for println, with predefined formatting.

# Input
- `fh::IO`:    (optional) file handle. If not supplied, prints to screen
- `z::String`: string, numbers and arrays to print

Paul.Soderlind@unisg.ch

"""
function printlnPs(fh::IO,z...;width=10,prec=3)

  for x in z                              #loop over inputs in z...
    if isa(x,AbstractArray)
      iob = IOBuffer()
      for i = 1:length(x)
        print(iob,fmtNumPs(x[i],width,prec,"right"))
      end
      print(fh,String(take!(iob)))
    else
      print(fh,fmtNumPs(x,width,prec,"right"))
    end
  end

  print(fh,"\n")

end
                      #when fh is not supplied: printing to screen
printlnPs(z...;width=10,prec=3) = printlnPs(stdout::IO,z...;width,prec)
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    fmtNumPs(z,width=10,prec=2,justify="right";prefix="",suffix="")

Create a formatted string of a float (eg, "%10.4f"), nothing (""),
while other values are passed through. Strings are right (or left) justified
and can optionally be given prefix and suffix (eg, ",")

# Notice
- With prec > 0 and isa(z,Integer), then the string is padded with 1+prec spaces
to align with the printing of floats with the same prec.

# Requires
- Printf (for 1.6-), fmtNumPsC (for < 1.6)

"""
function fmtNumPs(z,width=10,prec=2,justify="right";prefix="",suffix="")

  isa(z,Bool) && (z = convert(Int,z))             #Bool -> Int

  if isa(z,AbstractFloat)                         #example: 101.0234, prec=3
    if VERSION < v"1.6-"
      fmt    = "%$(width).$(prec)f"
      zRound = round(z,digits=prec)
      strLR  = fmtNumPsC(fmt,zRound)                #C fallback solution
    else
      fmt   = Printf.Format("%$(width).$(prec)f")
      strLR = Printf.format(fmt,z)
    end
  elseif isa(z,Nothing)
    strLR = ""
  elseif isa(z,Integer) && prec > 0               #integer followed by (1+prec spaces)
    strLR = string(z," "^(1+prec))
  else                                            #Int, String, Date, Missing, etc
    strLR = string(z)
  end

  strLR = string(prefix,strLR,suffix)

  if justify == "left"                            #justification
    strLR = rpad(strLR,width+length(prefix)+length(suffix))
  else
    strLR = lpad(strLR,width+length(prefix)+length(suffix))
  end

  return strLR

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    fmtNumPsC(fmt,z)

c fallback solution for formatting of floating point number. Used if VERSION < v"1.6-"
"""
function fmtNumPsC(fmt,z)                           #c fallback solution
  if ismissing(z) || isnan(z) || isinf(z)    #asprintf does not work for these cases
    str = string(z)
  else
    strp = Ref{Ptr{Cchar}}(0)
    len = ccall(:asprintf,Cint,(Ptr{Ptr{Cchar}},Cstring,Cdouble...),strp,fmt,z)
    str = unsafe_string(strp[],len)
    Libc.free(strp[])
  end
  return str
end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function printblue(x...)
  foreach(z->printstyled(z,color=:blue,bold=true),x)
  print("\n")
end
function printred(x...)
  foreach(z->printstyled(z,color=:red,bold=true),x)
  print("\n")
end
function printmagenta(x...)
  foreach(z->printstyled(z,color=:magenta,bold=true),x)
  print("\n")
end
function printyellow(x...)
  foreach(z->printstyled(z,color=:yellow,bold=true),x)
  print("\n")
end
#------------------------------------------------------------------------------

max_stars_repo_name: Japhiolite/JuliaTutorial

max_issues_repo_name: Japhiolite/JuliaTutorial

max_forks_repo_name: Japhiolite/JuliaTutorial


text:
using Pkg
Pkg.activate(pwd())

# # Wave equation
# ## Solving the wave equation

struct Wave
    f
    g
    c
end

# ### Exercise:
#
# Write the function `solve_wave(T, L, wave::Wave; n_t=100, n_x=100)` that solves the wave
# equation.
#
# **Hint**: Follow the procedure from the previous exercise. Discretize time and space,
# initialize the solution, add the boundary conditions, add the initial conditions and
# finally, iterate over time.
#
# ---
# ### Solution:

function solve_wave(T, L, wave::Wave; n_t=100, n_x=100)
    ts = range(0, T; length=n_t)
    xs = range(0, L; length=n_x)
    Δt = ts[2] - ts[1]
    Δx = xs[2] - xs[1]
    y = zeros(n_t, n_x)

    # boundary conditions
    y[:,1] .= wave.f(0)
    y[:,end] .= wave.f(L)

    # initial conditions
    y[1,2:end-1] = wave.f.(xs[2:end-1])
    y[2,2:end-1] = y[1,2:end-1] + Δt*wave.g.(xs[2:end-1])

    # solution for t = 2Δt, 3Δt, ..., T
    for t in 2:n_t-1, x in 2:n_x-1
        ∂y_xx = (y[t, x+1] - 2*y[t, x] + y[t, x-1])/Δx^2
        y[t+1, x] = c^2 * Δt^2 * ∂y_xx  + 2*y[t, x] - y[t-1, x]
    end

    return y
end

# ---

using Plots

function plot_wave(y, file_name; fps = 60, kwargs...)
    anim = @animate for (i, y_row) in enumerate(eachrow(y))
        plot(
            y_row;
            title = "t = $(i-1)Δt",
            xlabel = "x",
            ylabel = "y(t, x)",
            legend = false,
            linewidth = 2,
            kwargs...
        )
    end
    gif(anim, file_name; fps, show_msg = false)

    return nothing
end

# ### Exercise:
#
# Solve the wave equation for $L=\frac32\pi$, $T=240$, $c=0.02$ and the initial
# conditions
#
# $$
# \begin{aligned}
# f(x) &= 2e^{-(x-\frac L2)^2} + \frac{x}{L}, \\
# g(x) &= 0.
# \end{aligned}
# $$
#
# Use time discretization with stepsize $\Delta t=1$ and the space discretization with
# number of points $n_x=101$ and $n_x=7$ steps. Plot two graphs.
#
# ---
# ### Solution:

f(x,L) = 2*exp(-(x-L/2)^2) + x/L
g(x) = 0

L = 1.5*pi
T = 240
c = 0.02

wave = Wave(x -> f(x,L), g, c)

#+

y1 = solve_wave(T, L, wave; n_t=241, n_x=101)
plot_wave(y1, joinpath(pwd(), "lecture_13", "wave1.gif"); ylims=(-2,3), label="")

#+

y2 = solve_wave(T, L, wave; n_t=241, n_x=7)
plot_wave(y2, joinpath(pwd(), "lecture_13", "wave2.gif"); ylims=(-2,3), label="")


text:
##### Beginning of file

Test.@test( isdir(Snapshots.package_directory()) )

Test.@test( isdir(Snapshots.package_directory("ci")) )

Test.@test( isdir(Snapshots.package_directory("ci", "travis")) )

Test.@test( isdir(Snapshots.package_directory(TestModuleA)) )

Test.@test( isdir(Snapshots.package_directory(TestModuleB)) )

Test.@test(
    isdir( Snapshots.package_directory(TestModuleB, "directory2",) )
    )

Test.@test(
    isdir(
        Snapshots.package_directory(
            TestModuleB, "directory2", "directory3",
            )
        )
    )

Test.@test_throws(
    ErrorException,Snapshots.package_directory(TestModuleC),
    )

##### End of file


text:
module ConfigTest

using FactCheck

include(Pkg.dir("GenomeAnnotations", "src","Config.jl"))

conffile = Pkg.dir("GenomeAnnotations","conf","conf.ini")
if( isfile(conffile) )
   rm(conffile)
   rmdir( dirname(conffile ) )
end

LOCAL_STORE_PATH()

test_store_path = Pkg.dir("GenomeAnnotations","testdata","store", "annotation")
SET_LOCAL_STORE_PATH(test_store_path)

facts("local_store_path") do
    @fact LOCAL_STORE_PATH() => test_store_path
end

# set back to the original default
SET_LOCAL_STORE_PATH()

end


text:
"""
    root = Box{T}()
    child = Box(parent, level, edge)

Create a new box, a structure describing a rectangular domain nested
within a larger region. `parent` is the parent `Box`, `level` is an
integer assigned to indicate the likely priority that this box should be split,
and `edge` stores the value of one particular coordinate defining
the boundary between this box and its "partner" as defined by a golden
section split (or the edge of the outer bounds if this box lies between an
evaluation point and the outer bounds).

The clearest way to understand the design of `Box`
is by reference to Fig 2 in the [mcs](@ref) citation.
"""
mutable struct Box{T}
    level::Int
    edge::T  # the edge not corresponding to one of the parent's xvalues (splits that occur between dots in Fig 2)
    parent_cindex::Int # of its parent's children, which one is this?
    splitdim::Int # the dimension along which this box has been split (0 for leaf-boxes)
    parent::Box{T}
    xvalues::Vector{T} # the values of x_splitdim at which f is evaluated
    fvalues::Vector{T} # the corresponding values of f
    children::Vector{Box{T}}

    function Box{T}() where T
        # Create the root box
        box = new{T}(1, zero(T), typemax(Int)-1, 0)  # the root box will always be split along dimension 1
        box.parent = box
        return box
    end
    function Box{T}(parent::Box, level::Int, edge::T) where T
        # Create a new child and store it in the parent
        if !isdefined(parent, :children)
            parent.children = Box{T}[]
        end
        parent_cindex = length(parent.children) + 1
        box = new{T}(level, edge, parent_cindex, 0, parent)
        push!(parent.children, box)
        box
    end
end

Box(parent::Box{T}, splitdim, level, edge) where T =
    Box{T}(parent, splitdim, level, edge)

isroot(box::Box) = box.parent == box
issplit(box::Box) = isdefined(box, :children)


text:
export ode_order_lowering

function lower_varname(var::Variable, idv, order)
    order == 0 && return var
    name = Symbol(var.name, :_, string(idv.name)^order)
    return Variable(name; known = var.known)
end

function ode_order_lowering(sys::ODESystem)
    eqs_lowered, _ = ode_order_lowering(sys.eqs, sys.iv)
    ODESystem(eqs_lowered, sys.iv, [eq.x for eq in eqs_lowered], sys.ps)
end
function ode_order_lowering(eqs, iv)
    var_order = Dict{Variable,Int}()
    vars = Variable[]
    new_eqs = similar(eqs, DiffEq)
    new_vars = Variable[]

    for (i, eq) ∈ enumerate(eqs)
        var, maxorder = eq.x, eq.n
        if maxorder > get(var_order, var, 0)
            var_order[var] = maxorder
            any(isequal(var), vars) || push!(vars, var)
        end
        var′ = lower_varname(eq.x, iv, eq.n - 1)
        rhs′ = rename(eq.rhs)
        new_eqs[i] = DiffEq(var′, 1, rhs′)
    end

    for var ∈ vars
        order = var_order[var]
        for o in (order-1):-1:1
            lvar = lower_varname(var, iv, o-1)
            rvar = lower_varname(var, iv, o)
            push!(new_vars, rvar)

            rhs = rvar(iv())
            eq = DiffEq(lvar, 1, rhs)
            push!(new_eqs, eq)
        end
    end

    return (new_eqs, new_vars)
end

function rename(O::Expression)
    isa(O, Operation) || return O
    if is_derivative(O)
        (x, t, order) = flatten_differential(O)
        return lower_varname(x.op, t.op, order)(x.args...)
    end
    return Operation(O.op, rename.(O.args))
end


text:
F = Real;  V = Vector{<:F}

macro testrpositive()
    esc( :( if r == 0 && x == 0 throw(DomainError("x=r=0.0", "zero divided by zero")) end ) )
end

function ugly(x::F)
    1.0 / ( x^2 - 1.0 )
end

# scalar functions without arguments
"""
    bump(x::Real)

Compute a bump function at x.
"""
function bump(x::F)
    ifelse( abs(x) < 1.0, exp( ugly( x ) ), 0.0 )
end

"""
    dbump(x::Real)

Compute the first derivative of a bump function at x.
"""
function dbump(x::F)
    if abs(x) >= 1.0 return 0.0 end
    u = ugly(x)
    -2.0 * x * u^2 * exp(u)
end

"""
    ddbump(x::Real)

Compute the second derivative of a bump function at x.
"""
function ddbump(x::F)
    if abs(x) >= 1.0 return 0.0 end
    u = ugly(x)
    (6.0 * x^4 - 2.0) * u^4 * exp(u)
end

"""
    bump(x::Real, x0::Real, r::Real)

Compute a bump function at (x-x0) / r.
"""
function bump( x::F, x0::F, r::F )
    @testrpositive
    bump( ( x - x0 ) / r )
end

"""
    dbump(x::Real, x0::Real, r::Real)

Compute the first derivative of a bump function at (x-x0) / r.
"""
function dbump( x::F, x0::F, r::F )
    @testrpositive
    dbump( ( x - x0 ) / r ) / r
end

"""
    ddbump(x::Real, x0::Real, r::Real)

Compute the second derivative of a bump function at (x-x0) / r.
"""
function ddbump( x::F, x0::F, r::F )
    @testrpositive
    ddbump( ( x- x0 ) /r ) / r^2
end

# vector functions

function sumsqr( x::V )
    sum( x[i]^2 for i = 1: length(x) )
end

function ugly( x::V )
    1.0 / ( sumsqr( x ) - 1.0 )
end

"""
    bump(x::Vector{Real})

Compute a bump function at the vector x.
"""
function bump( x::V )
    ifelse( sumsqr(x) < 1.0, exp( ugly( x ) ), 0.0 )
end

"""
    dbump(x::Vector{Real})

Compute the first derivative of a bump function at the vector x.
"""
function dbump( x::V )
    if sumsqr(x) >= 1.0 return 0.0 end
    u = ugly( x )
    - 2.0 * x * u^2 * exp( u )
end

"""
    ddbump(x::Vector{Real})

Compute the second derivative of a bump function at the vector x.
"""
function ddbump( x::V )
    if sumsqr(x) >= 1.0 return 0.0 end
    u = ugly( x )
    [ 4.0 * u * ( 2.0 + u ) * x[i] * x[j] - ifelse( i == j, 2.0, 0.0 )
            for i = 1:length(x), j = 1:length(x) ] * u^2 * exp( u )
end

macro sanitycheck()
    esc(
        quote
            if length(x0) ≠ length(x) throw( DimensionMismatch("dimensions of x and x0 should match") ) end
            if length(r) ≠ length(x)  throw( DimensionMismatch("dimensions of x and r should match") ) end
            for i = 1:length(x)
                if r[i] == x[i] == 0.0 throw( DomainError("x=r=0.0", "zero divided by zero")) end
            end
        end )
end

"""
    bump(x::Vector{Real}, x0::Vector{Real}, r::Vector{Real})

Compute a vector bump function at x centered at x0 scaled by r (element-wise)
"""
function bump( x::V, x0::V, r::V )
    @sanitycheck
    bump( ( x - x0 ) ./ r )
end

"""
    dbump(x::Vector{Real}, x0::Vector{Real}, r::Vector{Real})

Compute the gradient of a vector bump function at x centered at x0 scaled by r (element-wise)
"""
function dbump( x::V, x0::V, r::V )
     @sanitycheck
    dbump( ( x - x0 ) ./ r ) ./ r
end

"""
    ddbump(x::Vector{Real}, x0::Vector{Real}, r::Vector{Real})

Compute the Hessian of a vector bump function at x centered at x0 scaled by r (element-wise)
"""
function ddbump( x::V, x0::V, r::V )
     @sanitycheck
    ddbump( ( x- x0 ) ./r ) ./ ( r .* r )
end


text:
#Electric field functions

function gaussian(amplitude, λ_central, τ_I_FWHM, τ, t)
    c = 299792458
    ω = 2 * π * c / λ_central
    E = @. amplitude * exp(-2*log(2)*((t - τ)/τ_I_FWHM)^2) * sin(ω*t)
    return E
end

function sinusoidal(amplitude, λ_central, t)
    c = 299792458
    ω = 2 * π * c / λ_central
    E = @. amplitude * sin(ω*t)
    return E
end


text:
include("run2.jl")
import Gadfly
Gadfly.set_default_plot_size(20Gadfly.cm, 12Gadfly.cm)

intrvls3 = Any[
    (0,2,50),
    (36.5,0.01,37)
]

# generate data for the plots on given intervals and save it to disk
function data3(intrvl)

    for i in 1:length(intrvl)
        qs = [intrvl[i][1]:intrvl[i][2]:intrvl[i][3]]
        Ps = zeros(length(qs))
        @sync @parallel for j in 1:length(qs)
            println("Finding ψ*ψ in bottom half for q = ",qs[j])
            Ps[j] = q3P(32,qs[j],20,1e-3);
        end
        writedlm("p3probs-$(intrvl[i][1])-to-$(intrvl[i][3]).dat",Ps)
    end

end

# create plots and save them as images
function plot3(intrvl)

    plots = Array(Gadfly.Plot,0)

    @sync @parallel for i in 1:length(intrvl)
        Ps = readdlm("p3probs-$(intrvl[i][1])-to-$(intrvl[i][3]).dat")
        p = Gadfly.plot(
            x=[intrvl[i][1]:intrvl[i][2]:intrvl[i][3]],
            y=Ps,
            Gadfly.Geom.line, Gadfly.Geom.point,
            Gadfly.Guide.xlabel("Charge of Point Particle, q"),
            Gadfly.Guide.ylabel("ψ*ψ"),
            Gadfly.Guide.title("Probability of Finding Particle in Bottom Half of Region vs. Charge of Point Particle"),
            Gadfly.Guide.xticks(ticks=[intrvl[i][1]:5intrvl[i][2]:intrvl[i][3]])
        )
        # save as SVG
        Gadfly.draw(Gadfly.SVG("p3probs-$(intrvl[i][1])-to-$(intrvl[i][3]).svg", 20Gadfly.cm, 12Gadfly.cm), p)
        push!(plots,p)

    end

    return plots

end


text:
using FactCheck
using FunctionalCollections

typealias PS PersistentSet

facts("Persistent Sets") do

    context("construction") do
        s = PS(1, 1, 2, 3, 3)
        @fact length(s) --> 3
        @fact length(PS{AbstractString}()) --> 0
        @fact typeof(PS{Integer}([1,2,3])) --> PS{Integer}
        @fact typeof(PS(Integer[1,2,3])) --> PS{Integer}
    end

    context("isequal") do
        @fact PS(1, 2, 3) --> PS(1, 2, 3)
        @fact PS(1, 2, 3) --> PS(3, 2, 1)
        @fact PS{AbstractString}() --> PS{Int}()
    end

    context("conj") do
        @fact conj(PS(1, 2, 3), 4) --> PS(1, 2, 3, 4)
        @fact conj(PS(1, 2, 3), 1) --> PS(1, 2, 3)
        @fact conj(PS(1, 2, 3), 4) --> PS(4, 3, 2, 1)
    end

    context("disj") do
        @fact disj(PS(1, 2, 3), 3) --> PS(1, 2)
        @fact disj(PS(1, 2), 3) --> PS(1, 2)
        @fact disj(PS{Int}(), 1234) --> PS{Int}()
    end

    context("in") do
        @fact "foo" in PS("foo", "bar") --> true
        @fact "baz" in PS("foo", "bar") --> false
    end

    context("filter") do
        @fact filter(iseven, PS(1, 2, 3, 4)) --> PS(2, 4)
    end

    context("setdiff, -") do
        @fact setdiff(PS(1, 2, 3), PS(1, 2)) --> PS(3)
        @fact setdiff(PS(1, 2), PS(1, 2, 3)) --> PS{Int}()
        @fact setdiff(PS(1, 2, 3), Set([1, 2])) --> PS(3)

        @fact PS(1, 2, 3) - PS(1, 2) --> PS(3)
    end

    context("isempty") do
        @fact PS{Int}() --> isempty
        @fact PS(1) --> not(isempty)
    end

    context("union") do
        @fact union(PS(1, 2, 3), PS(4, 5)) --> PS(1, 2, 3, 4, 5)
        @fact union(PS(1, 2, 3), PS(1, 2, 3, 4)) --> PS(1, 2, 3, 4)
        @fact union(PS(1), PS(2), PS(3)) --> PS(1, 2, 3)
    end

    context("isless, <=") do
        @fact PS(1)    <= PS(1, 2) --> true
        @fact PS(1, 2) <= PS(1, 2) --> true

        @fact PS(1, 2, 3) <= PS(1, 2) --> false
        @fact PS(1, 2) <= PS(1, 2, 3) --> true

        @fact isless(PS(1, 2), PS(1, 2))    --> false
        @fact isless(PS(1, 2), PS(1, 2, 3)) --> true
    end

    context("iteration") do
        @fact length([el for el in PS(1, 2, 3, 4)]) --> 4
    end

end


text:
importall MathProgBase.SolverInterface

mutable struct Problem

    alloc::Bool
    ptr::Ptr{Void}
    functions::Array{Function,1}
    constraints::Array{Constraint,1}
    net::Network

    function Problem(a, ptr, f, c, n)
        this = new(a, ptr, f, c, n)
        finalizer(this, dealloc)
        this
    end
end

Problem(net::Network) = Problem(true,
                                ccall((:PROB_new, libpfnet), Ptr{Void}, (Ptr{Void},), net.ptr),
                                Function[],
                                Constraint[],
                                net)

function add_constraint(prob::Problem, constr::Constraint)
    constr.alloc = false
    push!(prob.constraints, constr)
    ccall((:PROB_add_constr, libpfnet), Void, (Ptr{Void}, Ptr{Void},), prob.ptr, constr.ptr)
end

function add_function(prob::Problem, func::Function)
    func.alloc = false
    push!(prob.functions, func)
    ccall((:PROB_add_func, libpfnet), Void, (Ptr{Void}, Ptr{Void},), prob.ptr, func.ptr)
end

analyze(prob::Problem) = ccall((:PROB_analyze, libpfnet), Void, (Ptr{Void},), prob.ptr)

function apply_heuristics(prob::Problem, var_values::Base.Vector{Float64})
    vec::Ptr{Void} = pfnet_vec_from_Array(var_values)
    ccall((:PROB_apply_heuristics, libpfnet), Void, (Ptr{Void}, Ptr{Void},), prob.ptr, vec)
    Base.Libc.free(vec)
end

function combine_H(prob::Problem, coeff::Base.Vector{Float64})
    vec::Ptr{Void} = pfnet_vec_from_Array(coeff)
    ccall((:PROB_combine_H, libpfnet), Void, (Ptr{Void}, Ptr{Void}, Bool,), prob.ptr, vec, false)
    Base.Libc.free(vec)
end

function dealloc(prob::Problem)
    if prob.alloc
        ccall((:PROB_del, libpfnet), Void, (Ptr{Void},), prob.ptr)
    end
    prob.alloc = false
    prob.ptr = C_NULL
end

function eval(prob::Problem, var_values::Base.Vector{Float64})
    vec::Ptr{Void} = pfnet_vec_from_Array(var_values)
    ccall((:PROB_eval, libpfnet), Void, (Ptr{Void}, Ptr{Void},), prob.ptr, vec)
    Base.Libc.free(vec)
end

network(prob::Problem) = Network(ccall((:PROB_get_network, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr))

num_primal_variables(prob::Problem) = ccall((:PROB_get_num_primal_variables, libpfnet), Int, (Ptr{Void},), prob.ptr)
num_linear_equality_constraints(prob::Problem) = ccall((:PROB_get_num_linear_equality_constraints, libpfnet), Int, (Ptr{Void},), prob.ptr)
num_nonlinear_equality_constraints(prob::Problem) = ccall((:PROB_get_num_nonlinear_equality_constraints, libpfnet), Int, (Ptr{Void},), prob.ptr)

init_point(prob::Problem) = Vector(ccall((:PROB_get_init_point, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), true) # julia owns memory
upper_limits(prob::Problem) = Vector(ccall((:PROB_get_upper_limits, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), true) # julia owns memory
lower_limits(prob::Problem) = Vector(ccall((:PROB_get_lower_limits, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), true) # julia owns memory
x(prob::Problem) = Vector(ccall((:PROB_get_init_point, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), true) # julia owns memory

phi(prob::Problem) = ccall((:PROB_get_phi, libpfnet), Float64, (Ptr{Void},), prob.ptr)
gphi(prob::Problem) = Vector(ccall((:PROB_get_gphi, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), false)
Hphi(prob::Problem) = SparseMatrixCSC(ccall((:PROB_get_Hphi, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), false)

b(prob::Problem) = Vector(ccall((:PROB_get_b, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), false)
A(prob::Problem) = SparseMatrixCSC(ccall((:PROB_get_A, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), false)

f(prob::Problem) = Vector(ccall((:PROB_get_f, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), false)
J(prob::Problem) = SparseMatrixCSC(ccall((:PROB_get_J, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), false)

H_combined(prob::Problem) = SparseMatrixCSC(ccall((:PROB_get_H_combined, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr))

l(prob::Problem) = Vector(ccall((:PROB_get_l, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), false)
u(prob::Problem) = Vector(ccall((:PROB_get_u, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), false)
G(prob::Problem) = SparseMatrixCSC(ccall((:PROB_get_G, libpfnet), Ptr{Void}, (Ptr{Void},), prob.ptr), false)

show_problem(prob::Problem) = ccall((:PROB_show, libpfnet), Void, (Ptr{Void},), prob.ptr)

# min phi
# Ax = b
# f(x) = 0
# l <= Gx <= u
function solve(prob::Problem, solver::AbstractMathProgSolver)

    ma = size(A(prob))[1]
    mf = size(f(prob))[1]
    mg = size(G(prob))[1]
    nx = length(x(prob))
    ns = mg

    # min phi
    # 0 <= Ax-b <= 0
    # 0 <= f(x) - 0 <= 0
    # 0 <= Gx - s <= 0
    # [xmin l] <= [x s] <= [xmax u]

    numVar = nx + ns
    numConstr = ma+mf+mg
    l = [lower_limits(prob); pfnet.l(prob)]
    u = [upper_limits(prob); pfnet.u(prob)]
    lb = zeros(ma+mf+mg)
    ub = zeros(ma+mf+mg)
    sense = :Min
    model = NonlinearModel(solver)
    evaluator = ProblemEvaluator(prob, nx, ns, ma, mf, mg)
    y = [x(prob); zeros(ns)]

    loadproblem!(model, numVar, numConstr, l, u, lb, ub, sense, evaluator)
    setwarmstart!(model, y)
    optimize!(model)

    return (status(model),
            getsolution(model)[1:nx],
            getobjval(model))
end

mutable struct ProblemEvaluator <: AbstractNLPEvaluator
    problem::Problem
    nx::Int
    ns::Int
    ma::Int
    mf::Int
    mg::Int
end

function initialize(d::ProblemEvaluator, requested_features::Base.Vector{Symbol})
    analyze(d.problem)
end

function features_available(d::ProblemEvaluator)
    return Base.Vector([:Grad, :Jac, :Hess])
end

function eval_f(d::ProblemEvaluator, y)
    eval(d.problem, y[1:d.nx])
    return phi(d.problem)
end

function eval_g(d::ProblemEvaluator, g, y)
    x = y[1:d.nx]
    s = y[d.nx+1:end]
    eval(d.problem, x)
    g[:] = [A(d.problem)*x-b(d.problem);
            f(d.problem);
            G(d.problem)*x-s]
end

function eval_grad_f(d::ProblemEvaluator, g, y)
    eval(d.problem, y[1:d.nx])
    g[:] = [gphi(d.problem); zeros(d.ns)]
end

function jac_structure(d::ProblemEvaluator)
    Arow, Acol, Aval = SparseMatrixCOO(ccall((:PROB_get_A, libpfnet), Ptr{Void}, (Ptr{Void},), d.problem.ptr), false)
    Jrow, Jcol, Jval = SparseMatrixCOO(ccall((:PROB_get_J, libpfnet), Ptr{Void}, (Ptr{Void},), d.problem.ptr), false)
    Grow, Gcol, Gval = SparseMatrixCOO(ccall((:PROB_get_G, libpfnet), Ptr{Void}, (Ptr{Void},), d.problem.ptr), false)
    return ([Arow; Jrow+d.ma; Grow+d.ma+d.mf; (d.ma+d.mf+1):(d.ma+d.mf+d.mg)],
            [Acol; Jcol; Gcol; (d.nx+1):(d.nx+d.ns)])
end

function hesslag_structure(d::ProblemEvaluator)
    H1row, H1col, H1val = SparseMatrixCOO(ccall((:PROB_get_Hphi, libpfnet), Ptr{Void}, (Ptr{Void},), d.problem.ptr), false)
    H2row, H2col, H2val = SparseMatrixCOO(ccall((:PROB_get_H_combined, libpfnet), Ptr{Void}, (Ptr{Void},), d.problem.ptr), false)
    return([H1row; H2row],
           [H1col; H2col])
end

function eval_jac_g(d::ProblemEvaluator, J, y)
    x = y[1:d.nx]
    s = y[d.nx+1:end]
    eval(d.problem, x)
    Arow, Acol, Aval = SparseMatrixCOO(ccall((:PROB_get_A, libpfnet), Ptr{Void}, (Ptr{Void},), d.problem.ptr), false)
    Jrow, Jcol, Jval = SparseMatrixCOO(ccall((:PROB_get_J, libpfnet), Ptr{Void}, (Ptr{Void},), d.problem.ptr), false)
    Grow, Gcol, Gval = SparseMatrixCOO(ccall((:PROB_get_G, libpfnet), Ptr{Void}, (Ptr{Void},), d.problem.ptr), false)
    J[:] = [Aval; Jval; Gval; -ones(length(s))]
end

function eval_hesslag(d::ProblemEvaluator, H, y, sigma, mu)
    eval(d.problem, y[1:d.nx])
    combine_H(d.problem, mu[d.ma+1:d.ma+d.mf])
    H1row, H1col, H1val = SparseMatrixCOO(ccall((:PROB_get_Hphi, libpfnet), Ptr{Void}, (Ptr{Void},), d.problem.ptr), false)
    H2row, H2col, H2val = SparseMatrixCOO(ccall((:PROB_get_H_combined, libpfnet), Ptr{Void}, (Ptr{Void},), d.problem.ptr), false)
    H[:] = [sigma*H1val; H2val]
end


text:
const SPECIAL_PHS = Set{Symbol}([:x, :x1, :x2, :x3, :x4, :y, :z, :dz!dx, :dz!dy])

const SPECIAL_RULES = Dict{Any, Any}(
    # :(Main.foo) => [:(y = Main.foo(x)) => :(dz!dx = foo_grad(dz!dy, x))],

)

function special_derivative(g::EinGraph, iex::Expr, dzdy_vname::Symbol, z::Symbol, x::Symbol)
    ex = from_einstein(iex)
    op = ex.args[2].args[1]
    # in Einstein notation node may additionally contain sum(), unpack if so
    is_sum = false
    if op == :sum
        ex = :($(ex.args[1]) = $(ex.args[2].args[2]))
        op = ex.args[2].args[1]
        is_sum = true
    end
    rules = SPECIAL_RULES[op]
    dex = nothing
    for (pat, rpat) in rules
        wrt = split_deriv_name(rpat.args[1])[2] |> undname
        mt = matchex(pat, ex; phs=SPECIAL_PHS)
        if x == wrt && !isnull(mt)
            st = get(mt)
            st[:dz!dy] = dzdy_vname
            st[:dz!dx] = deriv_name(z, x)
            dex = subs(rpat, st)
            break
        end
    end
    # how to apply additional sum() here? should we include it into foo_grad!() itself?
    if dex != nothing
        lhs, rhs = dex.args
        @assert(getvalue(g[x]) != nothing, "$x is not evaluated in the source graph")
        lhs_w_idxs = with_indices(lhs, ndims(getvalue(g[x])))
        vnames = find_var_names(rhs)
        vars = [:($v[:]) for v in vnames]
        rhs_w_idxs = subs(rhs, Dict(vname => var for (vname, var) in zip(vnames, vars)))
        return TensorDeriv(g, :($lhs_w_idxs = $rhs_w_idxs))
    else
        error("Can't find special differentiation rule matching $iex")
    end
end

function _specialdiff(f_ex, df_ex)
    op = f_ex.args[2].args[1]
    dop = df_ex.args[2].args[1]
    push!(Espresso.SPECIAL_FUNCS, op, dop)
    if !haskey(SPECIAL_RULES, op)
        SPECIAL_RULES[op] = []
    end
    push!(SPECIAL_RULES[op], f_ex => df_ex)
end

macro specialdiff(f_ex, df_ex)
    _specialdiff(f_ex, df_ex)
    nothing
end


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

@testset "Quaternion.jl" begin

	 # function Quaternion()
   @testset "Quaternion" begin
      @test
      @test
      @test
      @test
   end

	 # function Quaternion(w::Float64,x::Float64,y::Float64,z::Float64)
   @testset "Quaternion" begin
      @test
      @test
      @test
      @test
   end

	 # function Quaternion(axis::Point3d,angle::Float64)
   @testset "Quaternion" begin
      @test
      @test
      @test
      @test
   end

	 # function convertToQuaternion(T::Matrix4)
   @testset "convertToQuaternion" begin
      @test
      @test
      @test
      @test
   end

	 # function convertToMatrix(q::Quaternion)
   @testset "convertToMatrix" begin
      @test
      @test
      @test
      @test
   end

end


text:
module Reversible

using Yao
using NiLang
using TropicalNumbers
using StaticArrays
using LinearAlgebra: Diagonal
using LuxurySparse
using ..TropicalYao: TropicalMatrixBlock, TropicalTypes, tropicalblock, copytensor, resettensor

export Tropical, CountingTropical, TropicalTypes
export LogLikeNumbers

include("LogLikeNumbers/LogLikeNumbers.jl")
using .LogLikeNumbers

include("instructs.jl")
include("spinglass_gates.jl")

end


text:
"""
module Fileaccess

This module contains low level fileaccess methods
"""
module Fileaccess

export load_jld_data,
       read_xyz,
       save_jld_data

using JLD

using ..Atoms
using ..Clusters
using ..Molecules

"""
save_jld_data(fname::AbstractString, data::Dict)

Saves given data to file in jld format.

Following keys are scanned from `data` and saved if present
- Method
- Basis
- cluster1
- cluster2
- Points
- Energy
- restart_energy
"""
function save_jld_data(fname::AbstractString, data::Dict)
    k = keys(data)
    jldopen(fname, "w") do file
        if haskey(data,"Method")
            file["Method"] = data["Method"]
            @info "Method information saved"
        end
        if haskey(data,"Basis")
            file["Basis"] = data["Basis"]
            @info "Basis informaiton saved"
        end
        if haskey(data,"cluster1")
            file["cluster1"] = data["cluster1"]
            @info "cluster1 information saved"
        end
        if haskey(data,"cluster2")
            file["cluster2"] = data["cluster2"]
            @info "cluster2 information saved"
        end
        if haskey(data,"Points")
            file["Points"] = data["Points"]
            @info "Points information saved"
        end
        if haskey(data,"Energy")
            file["Energy"] = data["Energy"]
            @info "Energy information saved"
        end
        if haskey(data,"restart_energy")
            file["restart_energy"] = data["restart_energy"]
            @info "Restart information saved"
        end
    end
    @info "Data writing to file \"$(fname)\" done"
end

"""
load_jld_data(fname::AbstractString)

Loads data from file and returns it as a [`Dict`](@ref)
"""
function load_jld_data(fname::AbstractString)
    data = load(fname)
    return data
end

"""
read_xyz(fname::AbstractString)

Reads xyz file and returns it as an Array of [`Cluster`](@ref)
"""
function read_xyz(fname::AbstractString)
    lines = Vector{String}()
    open(fname, "r") do file
        lines = readlines(file)
    end

    # How many atoms
    natoms = parse(Int, lines[1])

    # How many clusters - use of floor allows extra empty lines at end
    nclusters = Int(floor(length(lines)/(natoms+2)))
    @debug "Type of nclusters $(typeof(nclusters))"

    xyz = zeros(Float64, 3, natoms)
    atoms = Vector{AtomOnlySymbol}(undef, natoms)
    clusters = Vector{Cluster{AtomOnlySymbol}}()

    for nc in 1:nclusters
        for na in 1:natoms
            cont = split(lines[(nc-1)*(natoms+2)+na+2])
            atoms[na] = AtomOnlySymbol(cont[1])
            xyz[1,na] = parse(Float64, cont[2])
            xyz[2,na] = parse(Float64, cont[3])
            xyz[3,na] = parse(Float64, cont[4])
        end
        push!(clusters, deepcopy(Cluster(xyz, atoms)))
    end
    return clusters
end

end #module


text:
# Skipping MacroDefinition: IMPEXP __declspec ( dllimport )

const ngcomplex_t = Complex{Cdouble}

const VF_REAL = (1 << 0)
const VF_COMPLEX = (1 << 1)
const VF_ACCUM = (1 << 2)
const VF_PLOT = (1 << 3)
const VF_PRINT = (1 << 4)
const VF_MINGIVEN = (1 << 5)
const VF_MAXGIVEN = (1 << 6)
const VF_PERMANENT = (1 << 7)

"""
 Vector info obtained from any vector in ngspice.dll.
 Allows direct access to the ngspice internal vector structure,
 as defined in include/ngspice/devc.h .
"""
struct vector_info
    name::Cstring              # Same as so_vname
    type::Cint 	               # Same as so_vtype
    flags::Int16               # Flags (a combination of VF_*)
    realdata::Ptr{Cdouble}     # Real data
    compdata::Ptr{ngcomplex_t} # Complex data.
    length::Cint               # Length of the vector
end

const pvector_info = Ptr{vector_info}

struct vecvalues
    name::Cstring              # name of a specific vector
    creal::Cdouble             # actual data value
    cimag::Cdouble             # actual data value
    is_scale::Cint             # if 'name' is the scale vector
    is_complex::Cint           # if the data are complex numbers
end

const pvecvalues = Ptr{vecvalues}

struct vecvaluesall
    veccount::Cint         # number of vectors in plot
    vecindex::Cint         # index of actual set of vectors. i.e. the number of accepted data point
    vecsa::pvecvalues # values of actual set of vectors, indexed from 0 to veccount - 1
        #just pvecvalues?
end

const pvecvaluesall = Ptr{vecvaluesall}
struct vecinfo
    number::Cint           # number of vector, as postion in the linked list of vectors, s
    name::Cstring          # name of the actual vector
    is_real::Cint          # 1 if the actual vector has real data
    pdvec::Ptr{Cvoid}      # a void pointer to struct dvec *d, the actual vector
    pdvecscale::Ptr{Cvoid} # a void pointer to struct dvec *ds, the scale vector

end

const pvecinfo = Ptr{vecinfo}

struct vecinfoall
    name::Cstring
    title::Cstring
    date::Cstring
    type::Cstring
    veccount::Cint
    vecs::pvecinfo
end

#const pvecinfoall = Ptr{vecinfoall}

ngerrorf, bgrunningf, vecgetnum = 0, 0, 0
function std_print(text)
    occursin("stdout", text) &&
        (println(text[8:end]); return)
    (occursin("stderr", text) && !occursin("viewport", text)) &&
        (println(text[8:end]); return)
    occursin("viewport", text) &&
        (println("\nPlotting is skipped.\n")
                ; return)
    println(text)
end

function sendchar(_text::Ptr{Cchar}, id::Cint, userdata)::Cint
    _text != C_NULL || throw("Not a valid text")
    text = unsafe_string(_text)
    std_print(text)
    occursin(r"stderr Error:"i, text) && (ngerrorf = 1)
    return 0
end

gen_psendchar() = @cfunction(sendchar, Cint, (Ptr{Cchar}, Cint, Ptr{Cvoid}))

function sendstat(_text::Ptr{Cchar}, id::Cint, userdata)::Cint
    _text != C_NULL || throw("Not a valid text")
    text = unsafe_string(_text)
    std_print(text)
    return 0
end

gen_psendstat() = @cfunction(sendstat, Cint, (Ptr{Cchar}, Cint, Ptr{Cvoid}))

function bgthreadrunning(run::Cint, id::Cint, userdata::Ptr{Cvoid})::Cint
    bgrunning = run
    run ? println("BG thread is not running") :
        println("BG thread is running")
    return 0
end

gen_pbgthread() = @cfunction(bgthreadrunning, Cint, (Cint, Cint, Ptr{Cvoid}))

function controlledexit(exitstatus::Cint, immediate::Cint,
    quitexit::Cint, id::Cint, userdata::Ptr{Cvoid})::Cint
    quitexit == 1 && println("Returned from quit with exit status")
    immediate == 1 ? (println("Unloading NgSpice"); ngSpice_Command("quit")) :
        (println("Prepare an unload"); will_unload = 1)
    return exitstatus
end

gen_pcontrolledexit() = @cfunction(controlledexit, Cint, (Cint, Cint, Cint, Cint, Ptr{Cvoid}))

function senddata(vecdata::Ptr{vecinfoall},
    id::Cint, userdata::Ptr{Cvoid})::Cint
    return 0
end

gen_psenddata() = @cfunction(senddata, Cint, (Ptr{vecinfoall}, Cint, Ptr{Cvoid}))

function sendinitdata(initdata::Ptr{vecinfoall}, id::Cint, userdata::Ptr{Cvoid})
    #=
    This bit is problematic because it always returns C_NULLs even before
    `data.veccount` number of `vecinfo`s are passed.
    Even when that is handled it's o/p is always `nothing`
    and array of `nothing`s

    initdata == C_NULL && throw("No initialized data")
    data = unsafe_load(initdata)
    vec = unsafe_wrap(Array, data.vecs, (data.veccount, ))
    for v in vec
        v.name == C_NULL && return zero(Int32)
        vname = unsafe_string(v.name)
        vpdvec = unsafe_wrap(Array, v.pdvec, 10)
        vpdscale = unsafe_wrap(Array, v.pdvecscale, 10)
        println(vname, vpdscale, vpdvec)
    end
    =#
    return zero(Int32)
end

gen_psendinitdata() = @cfunction(sendinitdata, Cint, (Ptr{vecinfoall}, Cint, Ptr{Cvoid}))

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

_ntd_dist_and_shape(d::Distribution) = (d, varshape(d))

_ntd_dist_and_shape(s::ConstValueShape) = (ConstValueDist(s.value), s)

_ntd_dist_and_shape(s::IntervalSets.AbstractInterval) = _ntd_dist_and_shape(Uniform(minimum(s), maximum(s)))
_ntd_dist_and_shape(xs::AbstractVector{<:IntervalSets.AbstractInterval}) = _ntd_dist_and_shape(Product((s -> Uniform(minimum(s), maximum(s))).(xs)))
_ntd_dist_and_shape(xs::AbstractVector{<:Distribution}) = _ntd_dist_and_shape(Product(xs))
_ntd_dist_and_shape(x::Number) = _ntd_dist_and_shape(ConstValueShape(x))
_ntd_dist_and_shape(x::AbstractArray{<:Number}) = _ntd_dist_and_shape(ConstValueShape(x))

_ntd_dist_and_shape(h::Histogram{<:Real,1}) = _ntd_dist_and_shape(EmpiricalDistributions.UvBinnedDist(h))
_ntd_dist_and_shape(h::Histogram) = _ntd_dist_and_shape(EmpiricalDistributions.MvBinnedDist(h))

"""
    NamedTupleDist <: MultivariateDistribution
    NamedTupleDist <: MultivariateDistribution

A distribution with `NamedTuple`-typed variates.

`NamedTupleDist` provides an effective mechanism to specify the distribution
of each variable/parameter in a set of named variables/parameters.

Calling `varshape` on a `NamedTupleDist` will yield a
[`NamedTupleShape`](@ref).
"""
struct NamedTupleDist{
    names,
    DT <: (NTuple{N,Distribution} where N),
    AT <: (NTuple{N,ValueShapes.ValueAccessor} where N),
} <: Distribution{Multivariate,Continuous}
    _internal_distributions::NamedTuple{names,DT}
    _internal_shapes::NamedTupleShape{names,AT}
end

export NamedTupleDist

function NamedTupleDist(dists::NamedTuple{names}) where {names}
    dsb = map(_ntd_dist_and_shape, dists)
    NamedTupleDist(
        map(x -> x[1], dsb),
        NamedTupleShape(map(x -> x[2], dsb))
    )
end

@inline NamedTupleDist(;named_dists...) = NamedTupleDist(values(named_dists))

@inline _distributions(d::NamedTupleDist) = getfield(d, :_internal_distributions)
@inline _shapes(d::NamedTupleDist) = getfield(d, :_internal_shapes)

@inline Base.keys(d::NamedTupleDist) = keys(_distributions(d))

@inline Base.values(d::NamedTupleDist) = values(_distributions(d))

@inline function Base.getproperty(d::NamedTupleDist, s::Symbol)
    # Need to include internal fields of NamedTupleShape to make Zygote happy:
    if s == :_internal_distributions
        getfield(d, :_internal_distributions)
    elseif s == :_internal_shapes
        getfield(d, :_internal_shapes)
    else
        getproperty(_distributions(d), s)
    end
end

@inline function Base.propertynames(d::NamedTupleDist, private::Bool = false)
    names = propertynames(_distributions(d))
    if private
        (names..., :_internal_distributions, :_internal_shapes)
    else
        names
    end
end

ValueShapes.varshape(d::NamedTupleDist) = _shapes(d)

_ntd_length(d::Distribution) = length(d)
_ntd_length(d::ConstValueDist) = 0

Base.length(d::NamedTupleDist) = sum(_ntd_length, values(d))

function _ntd_logpdf(
    dist::ConstValueDist,
    acc::ValueShapes.ValueAccessor{<:ConstValueShape},
    x::AbstractVector{<:Real}
)
    float(zero(eltype(x)))
end

function _ntd_logpdf(
    dist::Distribution,
    acc::ValueShapes.ValueAccessor,
    x::AbstractVector{<:Real}
)
    logpdf(dist, float(x[acc]))
end

function Distributions.logpdf(d::NamedTupleDist, x::AbstractVector{<:Real})
    distributions = values(d)
    accessors = values(varshape(d))
    sum(map((dist, acc) -> _ntd_logpdf(dist, acc, x), distributions, accessors))
end

# ConstValueDist has no dof, so NamedTupleDist logpdf contribution must be zero:
_ntd_logpdf(dist::ConstValueDist, x::Any) = zero(Float32)

_ntd_logpdf(dist::Distribution, x::Any) = logpdf(dist, x)

function Distributions.logpdf(d::NamedTupleDist{names}, x::NamedTuple{names}) where names
    distributions = values(d)
    parvalues = values(x)
    sum(map((dist, d) -> _ntd_logpdf(dist, d), distributions, parvalues))
end

function _ntd_rand!(
    rng::AbstractRNG, dist::ConstValueDist,
    acc::ValueShapes.ValueAccessor{<:ConstValueShape},
    x::AbstractVector{<:Real}
)
    nothing
end

function _ntd_rand!(
    rng::AbstractRNG, dist::Distribution,
    acc::ValueShapes.ValueAccessor,
    x::AbstractVector{<:Real}
)
    rand!(rng, dist, view(x, acc))
    nothing
end

function Distributions._rand!(rng::AbstractRNG, d::NamedTupleDist, x::AbstractVector{<:Real})
    distributions = values(d)
    accessors = values(varshape(d))
    map((dist, acc) -> _ntd_rand!(rng, dist, acc, x), distributions, accessors)
    x
end

function _ntd_mode!(
    dist::ConstValueDist,
    acc::ValueShapes.ValueAccessor{<:ConstValueShape},
    params::AbstractVector{<:Real}
)
    nothing
end

function _ntd_mode!(
    dist::Distribution,
    acc::ValueShapes.ValueAccessor,
    params::AbstractVector{<:Real}
)
    view(params, acc) .= mode(dist)
    nothing
end

# Workaround, Distributions.jl doesn't define mode for Product:
function _ntd_mode!(
    dist::Distributions.Product,
    acc::ValueShapes.ValueAccessor,
    params::AbstractVector{<:Real}
)
    view(params, acc) .= map(mode, dist.v)
    nothing
end

function StatsBase.mode(d::NamedTupleDist)
    distributions = values(d)
    shape = varshape(d)
    accessors = values(shape)
    params = Vector{default_unshaped_eltype(shape)}(undef,shape)
    map((dist, acc) -> _ntd_mode!(dist, acc, params), distributions, accessors)
    params
end

function _ntd_var_or_cov!(A_cov::AbstractArray{<:Real,0}, dist::Distribution{Univariate})
    A_cov[] = var(dist)
    nothing
end

function _ntd_var_or_cov!(A_cov::AbstractArray{<:Real,2}, dist::Distribution{Multivariate})
    A_cov[:, :] = cov(dist)
    nothing
end

function _ntd_cov!(
    dist::ConstValueDist,
    acc::ValueShapes.ValueAccessor{<:ConstValueShape},
    A_cov::AbstractMatrix{<:Real}
)
    nothing
end

function _ntd_cov!(
    dist::Distribution,
    acc::ValueShapes.ValueAccessor,
    A_cov::AbstractMatrix{<:Real}
)
    _ntd_var_or_cov!(view(A_cov, acc, acc), dist)
    nothing
end

function Statistics.cov(d::NamedTupleDist)
    let n = length(d), A_cov = zeros(n, n)
        distributions = values(d)
        accessors = values(varshape(d))
        map((dist, acc) -> _ntd_cov!(dist, acc, A_cov), distributions, accessors)
        A_cov
    end
end


text:
using AtomsBase
using InteratomicPotentials
using StaticArrays
using Unitful
using UnitfulAtomic

B = zeros(length(0.1:0.1:2.5), 14)
dB = zeros(length(0.1:0.1:2.5), 14)
for (i, d) in enumerate(0.1:0.1:2.5)
    position1 = @SVector [0.0, 0.0, 0.0]
    position2 = @SVector [0.0, d, 0.0]

    element = :Ar
    atoms = [AtomsBase.Atom(element, position1 * 1u"Å"), AtomsBase.Atom(element, position2 * 1u"Å")]

    box = [[-4.0, 4.0], [-4.0, 4.0]]
    system = FlexibleSystem(box * 1u"Å", [Periodic(), Periodic()], atoms)

    snap = SNAPParams(5, 4, [:Ar], [4.0], 0.01, 0.989, [1.0])
    Btemp, dBtemp, W = compute_sna(system, snap)
    # print(Btemp)
    B[i, :] = sum(Btemp)
    dB[i, :] = dBtemp[1][:, 2]
end
show(stdout, "text/plain", B)
println(" ")

show(stdout, "text/plain", dB)

print(dB)

# show(stdout, "text/plain", dB)
# println(" ")

# show(stdout, "text/plain", W)
# println(" ")


text:
using Distributions
using Plots

# Sample x_{k+1} given x_k
# Parameters:
# - x_k: The kth X sample. One of -1, 0, or 1
# Outputs:
# - x_{k+1}: The (k+1)th X sample
function sample_x(xk)
    # Matrix defining p(x_{k+1}|x_k)
    probability_matrix = [0.9 0.05 0.05; 0.05 0.9 0.05; 0.05 0.05 0.9]

    # Extract row from probability matrix
    distribution = xk + 2

    # Define cumulative sum of distribution for sampling
    cumulative_distribution = cumsum(probability_matrix[Int(distribution),:])

    # Obtain x_k
    unif_sample = rand()
    xk_values = [-1 0 1]

    for i=1:3
        if unif_sample < cumulative_distribution[i]
            return xk_values[i]
        end
    end
end

# Obtains n samples of X based on HMM
# Inputs
# - n: An integer representing the number of samples to return
# Outputs
# A vector of length n where the kth element is the kth sample of X
function get_x_samples(n)
    # Store x samples
    samples = zeros(n)

    # Randomly set x1
    samples[1] = rand(-1:1)

    # Get samples x_2 through x_n
    for i=2:n
        samples[i] = sample_x(samples[i-1])
    end

    return samples
end

# Obtains n samples of Y based on HMM
# Inputs
# - x_samples: A vector of length n representing the samples of X
# Outputs
# A vector of length n where the kth element is the kth sample of Y
function get_y_samples(x_samples)
    # Store y samples
    n = length(x_samples)
    y_samples = zeros(n)

    for i=1:n
        𝒩 = Normal(x_samples[i],1)
        y_samples[i] = rand(𝒩)
    end

    return y_samples
end

# Obtian and plot samples
n = 1000
x_samples = get_x_samples(n)
y_samples = get_y_samples(x_samples)
plot(1:n,x_samples, seriestype=:scatter, bg = RGB(247/255, 236/255, 226/255), color = RGB(0,191/255,255/255), label = "x", alpha = 0.2)
plot!(1:n,y_samples, seriestype=:scatter, color = RGB(191/255,1,0), label =
"y", alpha = 0.2)

cd("/Users/elvis/Documents/DSI/courses/content/graphical-models/images/")
savefig("bn_sampling_ex_1.svg")

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
using IntUtils

function main()
    if length(ARGS) > 0
        n = parse(Int64, ARGS[1])
        a = Dict{Int64,Int64}[]
        for i in 2:n
            push!(a, prime_factors_table(i))
        end
        prod = 1
        for i in 2:n
            maxp = 0
            for j in 1:(n-1)
                p = get(a[j], i, 0)
                if p > maxp
                    maxp = p
                end
            end
            if maxp > 0
                prod *= i ^ maxp
            end
        end

        @printf("The smallest number that can be divided by each of the numbers from 1 to %d without any remainder is %d\n", n, prod)
    else
        error("this program needs a single numeric argument")
    end

end

main()


text:
using ArrayLayouts, BlockBandedMatrices, LinearAlgebra, Random, Test
import BlockBandedMatrices: colblockbandwidths,
    BroadcastStyle, BlockSkylineStyle, blocksizes

Random.seed!(0)

@testset "BlockSkylineMatrix" begin
    @testset "@jagot lmul! bug" begin
        rows = rand(1:10, 5)
        l = rand(-2:2, 5)
        u = rand(-2:2, 5)

        m = sum(rows)

        A = BlockSkylineMatrix(Zeros(m,m), (rows,rows), (l,u))
        @test @inferred(size(A)) == (sum(rows),sum(rows))

        @test Array(A) == Matrix(A)

        A.data .= rand(size(A.data)...)

        V = zeros(m,2)
        V[:,1] .= rand(m)
        reference = A*V[:,1]

        @view(V[:,2]) .= MulAdd(A, @view(V[:,1]))
        @test V[:,2] ≈ reference

        V[:,2] .= NaN
        @view(V[:,2]) .= MulAdd(A, @view(V[:,1]))
        @test V[:,2] ≈ reference
    end

    @testset "BlockSkylineMatrix multiplication" begin
        rows = [3, 1, 2, 1, 2, 1, 2, 1, 2, 1, 3]
        l,u = [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1], [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1]

        M = BlockSkylineMatrix{Float64}(undef, (rows,rows), (l,u))
        M.data .= 1

        d = Diagonal(1.0:size(M,2))
        D = BandedBlockBandedMatrix(d, (rows,rows), (0,0), (0,0))

        MD = M*D
        @test MD isa BlockSkylineMatrix
        @test MD == Matrix(M)*d
        @test colblockbandwidths(MD) == (l,u)

        MM = M*M
        @test MM isa BlockSkylineMatrix
        @test MM == Matrix(M)^2
        # Ensure correct (minimal) bandedness of product
        MMl,MMu = colblockbandwidths(MM)
        @test MMl[1:7] == [3,4,3,4,3,4,3]
        @test all(MMl[8:10] .≥ [3,2,1])
        @test all(MMu[2:4] .≥ [1,2,3])
        @test MMu[5:11] == [3,4,3,4,3,4,3]

        N = BlockBandedMatrix{Float64}(undef, (rows,rows), (1,1))
        N.data .= 1
        NN = N*N
        # We don't want a BlockBandedMatrix^2 to become a general
        # BlockSkylineMatrix
        @test NN isa BlockBandedMatrix
        @test NN == Matrix(N)^2

        rows = [9, 4, 1, 10, 6]
        O = BlockSkylineMatrix{Int64}(undef, (rows,rows), ([-2, 2, 0, 2, -1],[-1, 2, 1, 0, -1]))
        O.data .= 1
        OO = O*O
        @test OO isa BlockSkylineMatrix
        @test OO == Matrix(O)^2
    end

    @testset "Broadcasting" begin
        rows = [3, 1, 2, 1, 2, 1, 2, 1, 2, 1, 3]
        l,u = [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1], [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1]

        M = BlockSkylineMatrix{Float64}(undef, (rows,rows), (l,u))
        M.data .= 1

        @test BroadcastStyle(BlockSkylineMatrix) == BlockSkylineStyle()
        @test blocksizes(M + M) == blocksizes(M)

        # A is a BlockSkylineMatrix with varying bandwidths
        rows = [3,4,5]
        A = BlockSkylineMatrix{Float64}(undef, (rows,rows), ([1,0,1],[0,1,0]))
        A.data .= 1

        # B is diagonal
        B = BandedBlockBandedMatrix{Float64}(undef, (rows, rows), (0,0), (0,0))
        B.data .= -1

        # AB should have the same structure as A
        AB = A+B
        @test Matrix(AB) == Matrix(A) + Matrix(B)
        @test blocksizes(AB) == blocksizes(A)

        # C has larger bandwidths than A
        C = BandedBlockBandedMatrix{Float64}(undef, (rows, rows), (1,1), (1,1))
        C.data .= -1

        # Thus, AC should be a BlockSkylineMatrix with C's bandwidths.
        AC = A+C
        @test Matrix(AC) == Matrix(A) + Matrix(C)
        @test blocksizes(AC) == blocksizes(C)
    end
end