dharr

Dr. David Harrington

6636 Reputation

21 Badges

20 years, 80 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are replies submitted by dharr

@dharr Now the full VF2++ algorithm is implemented. Good for relatively sparse graphs but dense ones still difficult. Perhaps VF3...

Notes: It's currently a memory hog. Information about states that won't be visited could be deleted; T1 etc sets could be incrementally generated as needed. But for dense graphs, the VF3 algorithm may be better, since it (dynamically) uses information about G2 to prune the search tree; VF2++'s preprocessing isonly for G1.

Find if G1 is isomorphic to an induced subgraph of G2, using the VF2++ algorithm.
Graphs must be undirected without self-loops.

See VF2++ paper doi: 10.1016/j.dam.2018.02.018 for notation and description.
Code is in startup region (and is reproduced at the end of this worksheet)

restart;

with(GraphTheory):

claw:=CompleteGraph(1,3):
g:=CompleteGraph(2,2,2,2):
#DrawGraph(claw,size=[200,200]);
#DrawGraph(g,size=[200,200]);
IsSubgraphInducedIsomorphic(claw,g);

false

G1:=PathGraph(3):
G2:=Graph([1,2,3]):AddEdge(G2,{{3,1},{1,2}}):
DrawGraph(G1,size=[200,200]);
DrawGraph(G2,size=[200,200]);

IsSubgraphInducedIsomorphic(G1,G2,matching=true);

true, {1 = 2, 2 = 1, 3 = 3}

Snake in the box - see Wikipedia 

G1:=PathGraph(14); # longest for HyperGraph(5);
G2:=SpecialGraphs:-HypercubeGraph(5);
#DrawGraph(G1);
#DrawGraph(G2);

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14], Array(1..14, {(1) = {2}, (2) = {1, 3}, (3) = {2, 4}, (4) = {3, 5}, (5) = {4, 6}, (6) = {5, 7}, (7) = {6, 8}, (8) = {7, 9}, (9) = {8, 10}, (10) = {9, 11}, (11) = {10, 12}, (12) = {11, 13}, (13) = {12, 14}, (14) = {13}}), `GRAPHLN/table/15`, 0)

GRAPHLN(undirected, unweighted, ["00000", "00001", "00010", "00011", "00100", "00101", "00110", "00111", "01000", "01001", "01010", "01011", "01100", "01101", "01110", "01111", "10000", "10001", "10010", "10011", "10100", "10101", "10110", "10111", "11000", "11001", "11010", "11011", "11100", "11101", "11110", "11111"], Array(1..32, {(1) = {2, 3, 5, 9, 17}, (2) = {1, 4, 6, 10, 18}, (3) = {1, 4, 7, 11, 19}, (4) = {2, 3, 8, 12, 20}, (5) = {1, 6, 7, 13, 21}, (6) = {2, 5, 8, 14, 22}, (7) = {3, 5, 8, 15, 23}, (8) = {4, 6, 7, 16, 24}, (9) = {1, 10, 11, 13, 25}, (10) = {2, 9, 12, 14, 26}, (11) = {3, 9, 12, 15, 27}, (12) = {4, 10, 11, 16, 28}, (13) = {5, 9, 14, 15, 29}, (14) = {6, 10, 13, 16, 30}, (15) = {7, 11, 13, 16, 31}, (16) = {8, 12, 14, 15, 32}, (17) = {1, 18, 19, 21, 25}, (18) = {2, 17, 20, 22, 26}, (19) = {3, 17, 20, 23, 27}, (20) = {4, 18, 19, 24, 28}, (21) = {5, 17, 22, 23, 29}, (22) = {6, 18, 21, 24, 30}, (23) = {7, 19, 21, 24, 31}, (24) = {8, 20, 22, 23, 32}, (25) = {9, 17, 26, 27, 29}, (26) = {10, 18, 25, 28, 30}, (27) = {11, 19, 25, 28, 31}, (28) = {12, 20, 26, 27, 32}, (29) = {13, 21, 25, 30, 31}, (30) = {14, 22, 26, 29, 32}, (31) = {15, 23, 27, 29, 32}, (32) = {16, 24, 28, 30, 31}}), `GRAPHLN/table/16`, 0)

t,m:=IsSubgraphInducedIsomorphic(G1,G2,matching=true);
IsIsomorphic(G1,InducedSubgraph(G2,map(rhs,m)));

true, {1 = "00010", 2 = "00000", 3 = "00001", 4 = "00101", 5 = "01101", 6 = "01100", 7 = "01110", 8 = "11110", 9 = "10110", 10 = "10111", 11 = "10011", 12 = "11011", 13 = "11001", 14 = "11000"}

true

Coil in a box

G1:=CycleGraph(14); # longest for HyperGraph(5);
G2:=SpecialGraphs:-HypercubeGraph(5);
#DrawGraph(G1);
#DrawGraph(G2);

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14], Array(1..14, {(1) = {2, 14}, (2) = {1, 3}, (3) = {2, 4}, (4) = {3, 5}, (5) = {4, 6}, (6) = {5, 7}, (7) = {6, 8}, (8) = {7, 9}, (9) = {8, 10}, (10) = {9, 11}, (11) = {10, 12}, (12) = {11, 13}, (13) = {12, 14}, (14) = {1, 13}}), `GRAPHLN/table/21`, 0)

GRAPHLN(undirected, unweighted, ["00000", "00001", "00010", "00011", "00100", "00101", "00110", "00111", "01000", "01001", "01010", "01011", "01100", "01101", "01110", "01111", "10000", "10001", "10010", "10011", "10100", "10101", "10110", "10111", "11000", "11001", "11010", "11011", "11100", "11101", "11110", "11111"], Array(1..32, {(1) = {2, 3, 5, 9, 17}, (2) = {1, 4, 6, 10, 18}, (3) = {1, 4, 7, 11, 19}, (4) = {2, 3, 8, 12, 20}, (5) = {1, 6, 7, 13, 21}, (6) = {2, 5, 8, 14, 22}, (7) = {3, 5, 8, 15, 23}, (8) = {4, 6, 7, 16, 24}, (9) = {1, 10, 11, 13, 25}, (10) = {2, 9, 12, 14, 26}, (11) = {3, 9, 12, 15, 27}, (12) = {4, 10, 11, 16, 28}, (13) = {5, 9, 14, 15, 29}, (14) = {6, 10, 13, 16, 30}, (15) = {7, 11, 13, 16, 31}, (16) = {8, 12, 14, 15, 32}, (17) = {1, 18, 19, 21, 25}, (18) = {2, 17, 20, 22, 26}, (19) = {3, 17, 20, 23, 27}, (20) = {4, 18, 19, 24, 28}, (21) = {5, 17, 22, 23, 29}, (22) = {6, 18, 21, 24, 30}, (23) = {7, 19, 21, 24, 31}, (24) = {8, 20, 22, 23, 32}, (25) = {9, 17, 26, 27, 29}, (26) = {10, 18, 25, 28, 30}, (27) = {11, 19, 25, 28, 31}, (28) = {12, 20, 26, 27, 32}, (29) = {13, 21, 25, 30, 31}, (30) = {14, 22, 26, 29, 32}, (31) = {15, 23, 27, 29, 32}, (32) = {16, 24, 28, 30, 31}}), `GRAPHLN/table/22`, 0)

t,m:=IsSubgraphInducedIsomorphic(G1,G2,matching=true);
IsIsomorphic(G1,InducedSubgraph(G2,rhs~(m)));

true, {1 = "00000", 2 = "00001", 3 = "00101", 4 = "01101", 5 = "01100", 6 = "11100", 7 = "11000", 8 = "11001", 9 = "11011", 10 = "10011", 11 = "10111", 12 = "10110", 13 = "00110", 14 = "00010"}

true

But it is much slower to decide there are no matches when the cycle (or path) is one unit longer.

CodeTools:-Usage(IsSubgraphInducedIsomorphic(CycleGraph(15),G2,matching=true));

memory used=16.95GiB, alloc change=0.53GiB, cpu time=5.50m, real time=3.23m, gc time=2.75m

false

Random low density graph (length of output message here is a bug - SCR reported)

randomize(740472561204):
G2 := RandomGraphs:-RandomGraph(10000, 0.1);
NumberOfEdges(G2);
p:=combinat:-randcomb(Vertices(G2),30):
G1 := RelabelVertices(InducedSubgraph(G2,p),[$1..nops(p)]);

`[Length of output exceeds limit of 1000000]`

4999154

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], Array(1..30, {(1) = {}, (2) = {3, 23}, (3) = {2, 15, 17}, (4) = {7, 19, 20, 23, 30}, (5) = {12}, (6) = {8}, (7) = {4, 25, 26, 30}, (8) = {6, 25}, (9) = {17, 21, 23}, (10) = {17, 20}, (11) = {13, 21}, (12) = {5, 22, 25, 26}, (13) = {11, 24}, (14) = {19}, (15) = {3, 24, 25}, (16) = {20, 21, 27}, (17) = {3, 9, 10}, (18) = {21, 22, 23}, (19) = {4, 14, 29}, (20) = {4, 10, 16}, (21) = {9, 11, 16, 18, 28, 30}, (22) = {12, 18, 26}, (23) = {2, 4, 9, 18}, (24) = {13, 15, 25, 26}, (25) = {7, 8, 12, 15, 24}, (26) = {7, 12, 22, 24, 28}, (27) = {16, 29}, (28) = {21, 26}, (29) = {19, 27}, (30) = {4, 7, 21}}), `GRAPHLN/table/29`, 0)

t,m:=CodeTools:-Usage(IsSubgraphInducedIsomorphic(G1,G2,matching=true));
IsIsomorphic(G1,InducedSubgraph(G2,rhs~(m)));

memory used=483.58MiB, alloc change=-114.43MiB, cpu time=3.09s, real time=2.59s, gc time=718.75ms

true, {1 = 21, 2 = 1869, 3 = 1209, 4 = 139, 5 = 264, 6 = 34, 7 = 269, 8 = 16, 9 = 5, 10 = 844, 11 = 42, 12 = 28, 13 = 202, 14 = 167, 15 = 346, 16 = 13, 17 = 50, 18 = 20, 19 = 140, 20 = 55, 21 = 1, 22 = 197, 23 = 334, 24 = 190, 25 = 385, 26 = 272, 27 = 144, 28 = 74, 29 = 187, 30 = 29}

true

Random smaller higher density graph is much harder

randomize(22807444796519):
G2 := RandomGraphs:-RandomGraph(50, 0.5);
p:=combinat:-randcomb(Vertices(G2),11):
G1 := RelabelVertices(InducedSubgraph(G2,p),[$1..nops(p)]);

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50], Array(1..50, {(1) = {2, 4, 5, 8, 9, 10, 11, 14, 15, 16, 17, 20, 22, 24, 27, 28, 29, 30, 32, 33, 34, 37, 38, 41, 42, 43, 44, 45, 46, 47, 48}, (2) = {1, 7, 9, 10, 11, 14, 16, 17, 19, 20, 21, 27, 28, 31, 33, 39, 40, 41, 45, 47}, (3) = {4, 6, 7, 8, 9, 10, 12, 16, 17, 18, 24, 28, 29, 33, 34, 37, 38, 39, 40, 42, 44, 46, 47}, (4) = {1, 3, 5, 6, 8, 9, 10, 12, 18, 20, 21, 23, 24, 25, 27, 28, 29, 30, 34, 37, 38, 39, 40, 41, 44, 45, 46, 47, 48, 50}, (5) = {1, 4, 6, 8, 11, 12, 14, 15, 16, 17, 18, 22, 23, 24, 27, 29, 30, 32, 33, 34, 35, 36, 37, 38, 40, 42, 43, 49, 50}, (6) = {3, 4, 5, 7, 9, 10, 12, 15, 16, 18, 19, 21, 22, 23, 25, 27, 28, 29, 32, 34, 36, 37, 38, 40, 41, 42, 43, 44, 47, 50}, (7) = {2, 3, 6, 9, 10, 11, 13, 14, 16, 17, 18, 19, 20, 23, 24, 26, 29, 30, 31, 33, 37, 39, 42, 43, 45, 47}, (8) = {1, 3, 4, 5, 10, 11, 12, 13, 17, 20, 26, 27, 28, 30, 33, 34, 35, 36, 40, 42, 45, 46, 49, 50}, (9) = {1, 2, 3, 4, 6, 7, 17, 20, 21, 24, 25, 28, 29, 33, 36, 37, 38, 40, 43, 45, 49, 50}, (10) = {1, 2, 3, 4, 6, 7, 8, 11, 13, 14, 15, 18, 21, 23, 24, 26, 28, 30, 31, 34, 36, 38, 39, 40, 42, 44, 45}, (11) = {1, 2, 5, 7, 8, 10, 12, 13, 14, 16, 17, 19, 20, 23, 25, 26, 27, 28, 30, 31, 32, 34, 36, 44, 45, 46, 47, 50}, (12) = {3, 4, 5, 6, 8, 11, 15, 17, 20, 21, 22, 24, 25, 27, 28, 29, 30, 31, 34, 36, 37, 39, 42, 43, 44, 45, 46, 48, 49, 50}, (13) = {7, 8, 10, 11, 15, 17, 18, 19, 21, 23, 24, 26, 27, 29, 30, 33, 34, 35, 36, 39, 40, 41, 48}, (14) = {1, 2, 5, 7, 10, 11, 15, 18, 19, 20, 23, 24, 25, 29, 32, 34, 35, 36, 39, 40, 41, 45, 46, 47, 48, 49}, (15) = {1, 5, 6, 10, 12, 13, 14, 19, 22, 23, 24, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 41, 45, 47, 50}, (16) = {1, 2, 3, 5, 6, 7, 11, 20, 22, 24, 25, 26, 28, 32, 33, 34, 37, 42, 45, 47, 49}, (17) = {1, 2, 3, 5, 7, 8, 9, 11, 12, 13, 19, 21, 23, 24, 25, 26, 27, 28, 29, 30, 36, 37, 39, 42, 44, 45, 47, 49, 50}, (18) = {3, 4, 5, 6, 7, 10, 13, 14, 21, 23, 24, 26, 29, 30, 31, 33, 35, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 50}, (19) = {2, 6, 7, 11, 13, 14, 15, 17, 21, 22, 23, 25, 26, 27, 29, 30, 31, 32, 33, 36, 37, 38, 39, 41, 42, 43, 44, 45, 46, 49, 50}, (20) = {1, 2, 4, 7, 8, 9, 11, 12, 14, 16, 21, 22, 23, 25, 28, 29, 30, 31, 33, 35, 36, 39, 40, 41, 44, 45, 46, 47}, (21) = {2, 4, 6, 9, 10, 12, 13, 17, 18, 19, 20, 22, 23, 25, 27, 28, 31, 33, 37, 41, 43, 44, 45, 48, 49}, (22) = {1, 5, 6, 12, 15, 16, 19, 20, 21, 30, 33, 34, 36, 37, 39, 40, 41, 42, 45, 46, 48, 49, 50}, (23) = {4, 5, 6, 7, 10, 11, 13, 14, 15, 17, 18, 19, 20, 21, 24, 25, 26, 28, 29, 34, 39, 41, 42, 43, 45, 47, 48, 50}, (24) = {1, 3, 4, 5, 7, 9, 10, 12, 13, 14, 15, 16, 17, 18, 23, 27, 28, 32, 33, 37, 38, 40, 41, 44, 46, 49}, (25) = {4, 6, 9, 11, 12, 14, 16, 17, 19, 20, 21, 23, 26, 27, 29, 30, 33, 35, 37, 41, 45, 46, 49}, (26) = {7, 8, 10, 11, 13, 16, 17, 18, 19, 23, 25, 27, 29, 30, 34, 35, 36, 38, 39, 40, 42, 43, 45}, (27) = {1, 2, 4, 5, 6, 8, 11, 12, 13, 17, 19, 21, 24, 25, 26, 28, 29, 31, 32, 34, 35, 36, 38, 40, 41, 42, 45, 46, 47}, (28) = {1, 2, 3, 4, 6, 8, 9, 10, 11, 12, 15, 16, 17, 20, 21, 23, 24, 27, 30, 32, 33, 36, 37, 39, 40, 44, 46, 47, 49}, (29) = {1, 3, 4, 5, 6, 7, 9, 12, 13, 14, 15, 17, 18, 19, 20, 23, 25, 26, 27, 34, 43, 45, 48, 49, 50}, (30) = {1, 4, 5, 7, 8, 10, 11, 12, 13, 15, 17, 18, 19, 20, 22, 25, 26, 28, 32, 33, 34, 36, 39, 41, 42, 43, 46}, (31) = {2, 7, 10, 11, 12, 15, 18, 19, 20, 21, 27, 36, 38, 39, 40, 42, 44, 48, 49, 50}, (32) = {1, 5, 6, 11, 14, 15, 16, 19, 24, 27, 28, 30, 33, 36, 38, 39, 41, 42, 46, 47, 48, 50}, (33) = {1, 2, 3, 5, 7, 8, 9, 13, 16, 18, 19, 20, 21, 22, 24, 25, 28, 30, 32, 34, 35, 37, 38, 39, 40, 42, 44, 45, 50}, (34) = {1, 3, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 22, 23, 26, 27, 29, 30, 33, 37, 38, 39, 42, 44, 48, 49}, (35) = {5, 8, 13, 14, 15, 18, 20, 25, 26, 27, 33, 40, 41, 42, 43, 45, 48, 49, 50}, (36) = {5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 17, 19, 20, 22, 26, 27, 28, 30, 31, 32, 37, 40, 41, 47, 48, 49}, (37) = {1, 3, 4, 5, 6, 7, 9, 12, 15, 16, 17, 19, 21, 22, 24, 25, 28, 33, 34, 36, 42, 44, 45, 47, 48}, (38) = {1, 3, 4, 5, 6, 9, 10, 15, 18, 19, 24, 26, 27, 31, 32, 33, 34, 41, 43, 46, 48, 49}, (39) = {2, 3, 4, 7, 10, 12, 13, 14, 15, 17, 18, 19, 20, 22, 23, 26, 28, 30, 31, 32, 33, 34, 40, 41, 43, 45, 46, 48, 49}, (40) = {2, 3, 4, 5, 6, 8, 9, 10, 13, 14, 18, 20, 22, 24, 26, 27, 28, 31, 33, 35, 36, 39, 41, 42, 45, 48, 49, 50}, (41) = {1, 2, 4, 6, 13, 14, 15, 18, 19, 20, 21, 22, 23, 24, 25, 27, 30, 32, 35, 36, 38, 39, 40, 45, 46, 47, 48, 49, 50}, (42) = {1, 3, 5, 6, 7, 8, 10, 12, 16, 17, 18, 19, 22, 23, 26, 27, 30, 31, 32, 33, 34, 35, 37, 40, 44, 45, 48, 49, 50}, (43) = {1, 5, 6, 7, 9, 12, 18, 19, 21, 23, 26, 29, 30, 35, 38, 39, 49, 50}, (44) = {1, 3, 4, 6, 10, 11, 12, 17, 18, 19, 20, 21, 24, 28, 31, 33, 34, 37, 42, 45, 50}, (45) = {1, 2, 4, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 25, 26, 27, 29, 33, 35, 37, 39, 40, 41, 42, 44, 47, 48}, (46) = {1, 3, 4, 8, 11, 12, 14, 18, 19, 20, 22, 24, 25, 27, 28, 30, 32, 38, 39, 41, 47, 48, 49}, (47) = {1, 2, 3, 4, 6, 7, 11, 14, 15, 16, 17, 18, 20, 23, 27, 28, 32, 36, 37, 41, 45, 46, 50}, (48) = {1, 4, 12, 13, 14, 18, 21, 22, 23, 29, 31, 32, 34, 35, 36, 37, 38, 39, 40, 41, 42, 45, 46, 49}, (49) = {5, 8, 9, 12, 14, 16, 17, 19, 21, 22, 24, 25, 28, 29, 31, 34, 35, 36, 38, 39, 40, 41, 42, 43, 46, 48, 50}, (50) = {4, 5, 6, 8, 9, 11, 12, 15, 17, 18, 19, 22, 23, 29, 31, 32, 33, 35, 40, 41, 42, 43, 44, 47, 49}}), `GRAPHLN/table/34`, 0)

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], Array(1..11, {(1) = {3, 5, 7, 11}, (2) = {3, 4, 9}, (3) = {1, 2, 9, 10, 11}, (4) = {2, 5, 6, 7, 9, 11}, (5) = {1, 4, 9, 10, 11}, (6) = {4, 7, 9}, (7) = {1, 4, 6, 8, 9}, (8) = {7}, (9) = {2, 3, 4, 5, 6, 7, 10}, (10) = {3, 5, 9}, (11) = {1, 3, 4, 5}}), `GRAPHLN/table/35`, 0)

t,m:=CodeTools:-Usage(IsSubgraphInducedIsomorphic(G1,G2,matching=true));
IsIsomorphic(G1,InducedSubgraph(G2,rhs~(m)));

memory used=5.47GiB, alloc change=0 bytes, cpu time=54.72s, real time=47.40s, gc time=8.72s

true, {1 = 13, 2 = 46, 3 = 8, 4 = 25, 5 = 29, 6 = 37, 7 = 21, 8 = 2, 9 = 4, 10 = 50, 11 = 26}

true

NULL

Code from startup region reproduced here

IsSubgraphInducedIsomorphic := proc(G1::Graph, G2::Graph, {matching::truefalse:=false})
        description "Determines if undirected loopless G1 is isomorphic to an induced subgraph of G2",
                        "and if so optionally returns the matching function for the vertices";
        local n1,n2,D1,W1,V1,A1orig,A1,D2,W2,V2,A2,perm,MM,Morder,r,Vd,V1remaining,
                depth,allV1,allV2,Mout,Mstate,expandmatching,Premaining,
                M,M1,M2,T1,T1sorted,T2,T1t,T2t,P,pair,consistent,cut,i,j;
        D1,W1,V1,A1orig := op(1..4,G1); # directedness,weightedness,vertices,neighbors
        D2,W2,V2,A2     := op(1..4,G2);
        if D1 <> 'undirected' or D2 <> 'undirected' then error "both graphs must be undirected" end if;
        n1 := nops(V1):                # number of vertices in G1
        n2 := nops(V2):
        if orseq(i in A1orig[i], i = 1..n1) or orseq(i in A2[i], i=1..n2) then error "graph(s) have self-loops" end if;
        allV1 := {$1..n1};                # all vertices in G1
        allV2 := {$1..n2};
        # simple rejections based on number of vertices and max degrees  
        if n1 > n2 then return false end if;
        if max(map(nops,A1orig)) > max(map(nops,A2)) then return false end if;
        # First presort vertices of G1 in special BFS order
        # We are not using labels so can presort in decreasing degree order outside outermost while loop
        V1remaining := sort([$1..n1], (a,b) -> nops(A1orig[a]) >= nops(A1orig[b]));
        Morder:=Array([]);
        MM:={};
        while numelems(V1remaining) <> 0 do # loop for each connected component
                r := V1remaining[1]; # root for BFS tree for this component
                Vd := [r]; # depth = 0
                while nops(Vd) > 0 do # for each depth
                        # sort according to decreasing connections to M and
                        # otherwise decreasing degree
                        if nops(Vd) > 1 then
                                Vd:=sort(Vd, proc(a,b)
                                                        local na := nops(A1orig[a] intersect MM),
                                                                 nb := nops(A1orig[b] intersect MM);
                                                        na > nb  or na = nb and nops(A1orig[a]) >= nops(A1orig[b])
                                                   end proc)
                        end if;
                        # add these to the order and update
                        Morder ,= Vd[];
                        MM := MM union {Vd[]};
                        V1remaining := remove(has, V1remaining, Vd);
                        # go one level deeper by finding new neighbors of Vd;
                        Vd:=convert(`union`(seq(A1orig[Vd])) minus MM, list);
                end do;
        end do;        
        perm:=[seq(Morder)]; # order of vertices we want.
        # We now reorder and relabel the graph so the chosen order becomes the natural order
        A1 := A1orig[perm];
        A1 := eval(A1, perm =~ [$1..n1]);
        # Start of main part of algorithm
        # current matching M is represented as a set of equations
        # that map G1 vertices (lhs) to G2 vertices (rhs), e.g., (2=4,3=5}
        #
        # expandmatching works out M1,M2,T1,T1sorted,T2,T1t,T2t for a given M
        expandmatching := proc(M) option remember;
                local
                        M1 := map(lhs,M), # set of verts in G1 in matching
                        M2 := map(rhs,M), # set of verts in G2 in matching
                        T1 := {seq(op(A1[i]), i in M1)} minus M1, # set of neighbors of M1 not in M1
                        # sort T1 in order of decreasing numbers of neighbors in M1
                        T1sorted := sort([T1[]],(a,b)->nops(A1[a] intersect M1) >= nops(A1[b] intersect M1)),
                        T2 := {seq(op(A2[i]), i in M2)} minus M2, # set of neighbors of M2 not in M2
                        # T1t = T1tilde are further away vertices
                        T1t := allV1 minus M1 minus T1,
                        T2t := allV2 minus M2 minus T2;
                M1,M2,T1,T1sorted,T2,T1t,T2t
        end proc;
        # Premaining is table indexed by a matching that gives
        # the remaining addition possibilities (P) not yet tried for this matching
        Premaining := table('sparse');
        M := {};
        Mout := 0;
        depth := 0;
        # main loop - each iteration assesses a candidate state M
        # in a depth-first search
        while depth >= 0 do
                M1,M2,T1,T1sorted,T2,T1t,T2t := expandmatching(M);
                if M1 = allV1 then # success
                        Mout := M;
                        break
                end if;
                # find remaining candidate pairs for addition to this state
                if Premaining[M] = 0 then # new M
                        if nops(T1) <> 0 then
                                P := {seq([T1sorted[1],j],j in T2)};
                        elif         nops(allV1 minus M1) > 0 then        
                                P := {seq([(allV1 minus M1)[1],j], j in (allV2 minus M2))};
                        else
                                P:={}
                        end if;
                        Premaining[M] := P;
                else
                        P := Premaining[M];
                end if;
                # check out possibilities and go deeper if one is found
                while numelems(P) > 0 do
                        pair := P[1];
                        P := P minus {pair};
                        Premaining[M] := P;
                        i,j := pair[];
                        # consistency check for [i,j]
                        # neighbors of i in M1 must correspond to neighbors of j in M2
                        consistent:=evalb(eval(A1[i] intersect M1, M) = (A2[j] intersect M2));
                        # cut rule
                        cut := evalb(nops(A2[j] intersect T2) < nops(A1[i] intersect T1)
                                or nops(A2[j] intersect T2t) < nops(A1[i] intersect T1t));
                        if consistent and not cut then
                                Mstate[depth] := M; # save state
                                depth := depth + 1;  
                                M := M union {i = j};
                                next 2
                        end if
                end do;
                # no possibilities, backtrack
                depth := depth - 1;
                M := Mstate[depth];
        end do;
        forget(expandmatching);
        # return values
        if Mout <> 0 then   
                if matching then
                        true, map(eq->V1[perm[lhs(eq)]]=V2[rhs(eq)], Mout)
                else
                        true
                    end if;
          else
                false
        end if;
end proc:

NULL

Download VF2pp.mw

@mmcdara @jalal The "this:///..." notation refers to a file within the workbook (not a .mw file), where "this" means within the current workbook; see

this help

So I think the original problem was, as @mmcdara pointed out, that it was not a string. But you don't wand to append this to currentdir(). I think the Read needs only the "this///...", but until you upload your workbook it is hard for someone to check this. As @acer pointed out, you need to compress your workbook to a .zip file before it can be uploaded to Mapleprimes - I see you have now done that.

@sursumCorda Yes, it is very fast. Adding iterations=100 to the CodeTools:-Usage shows they are both about the same time, so the apparently much faster time is a fluctuation. Vectors of integers also can be Trimmed. 

@jalal 

NULL

NULL

restart

randomize()

A := RandomTools:-Generate(integer(range = -5 .. 5, exclude = {0}))

-3

B := RandomTools:-Generate(integer(range = -5 .. 5, exclude = {0, A}))

-2

operation := combinat:-randcomb([`+`, `-`, `*`], 1)[]

`*`

operation(A, B)

6

 

 

 

Download Qop.mw

@mmcdara I agree with the OPs observations of errors in Maple 2023. I found them hard to track down, and got confused by the gammas. But please check, you have both gamma_i and gamma__i in the gap procedure. .I usually define the ode and ic set outside of the procedure to be passed to minimize, and do eval on them within the procedure.

@mmcdara Your 

solution := tanh((1/4)*sqrt(2)*x)

is I believe the solution the OP referred to. It is an exact solution for z(0)=0, z(infinity)=1 (as verified by taking the limit), but is only an approximate solution for z(0)=0, z(15)=1.

@mmcdara It is a confusing situation - while there isn't a help page for, say, `[]`, asking for help does take you to a page that descibes the user-level equivalent. And yet there is a help page for `?[]`.

@mmcdara I think the OP wanted documentation for their usage in the form `<,>`(1,2,3) and `<|>`(1,2,3).

@mmcdara Nice analysis, vote up. While having the procedure name integraal the same as the local integraal may be confusing, it is not an error, as the procedure name is global and does not conflict with the local variable within the procedure.

ps: the loop within optimalisatie looks to me like

for i while i < 10 and 1/10 < abs(int__sim - int__ex) do

(see above)

I recently looked for the documentation for `<|>` and couldn't find any. For `<,>`( sequence ) the user level <sequence> is equivalent and turns a sequence into a column vector (or perhaps matrix if the sequence has more complicated entries). Likewise `<|>`(sequence) turns the sequence into a row vector as though you entered <a|b|c...>. Programmatically, one may have generated a sequence and then the user level <sequence> is fine, which may be the reson that `<,>` is undocumented. On the other hand `<|>` is the only way to turn a sequence into a row vector, so there is a better case for documentation here.

Of course I may have just failed to find the documentation also...

If I paste the optimalisatie code into a code edit region I get the message that maplemint got about "These local variables were assigned a value, but otherwise unused:   int__sim_min", which is correct: You calculate this quantity (which looks like int with subscript sim_min) but then do not use it.

In the code edit region there is no message "invalid left hand side of assignment", so that seems to be spurious. [Edit; as @mmcdara points out this is assigning to Gamma__i_min and Gamma__i_max]

Note that the printed code for the procedure has many things such as int__sim_min looking like int[sim_min] as though you entered them as indexed variables. But pasting the 2-D code into a 1-D region or converting it to 1-D with the context menu shows that actually they are single variables displaying with subscripts, as you presumably intended. I've never seen that before.

In a code edit region it looks like

optimalisatie := proc(t__ex, c__ex, C__0, Q, V__r, m__b, rho, R, Gamma__i_min, Gamma__i_max)
	local Gamma__i, t__sim, C__i, int__sim, int__ex, C__i_min, C__i_max, C__i0_min, C__i0_max, C__i0, int__sim_min, int__sim_max, i;
	Gamma__i := 1/2*Gamma__i_min + 1/2*Gamma__i_max;
	t__sim, C__i := CC(C__0, Q, V__r, m__b, rho, R, Gamma__i);
	int__sim := integraal(t__sim, C__i);
	int__ex := integraal(t__ex, c__ex);
	for i while i < 10 and 1/10 < abs(int__sim - int__ex) do
		Gamma__i := 1/2*Gamma__i_min + 1/2*Gamma__i_max;
		t__sim, C__i_min := CC(C__0, Q, V__r, m__b, rho, R, Gamma__i_min);
		t__sim, C__i_max := CC(C__0, Q, V__r, m__b, rho, R, Gamma__i_max);
		t__sim, C__i := CC(C__0, Q, V__r, m__b, rho, R, Gamma__i);
		C__i0_min := normalisatie(C__i_min, C__0);
		C__i0_max := normalisatie(C__i_max, C__0);
		C__i0 := normalisatie(C__i, C__0);
		int__sim_min := integraal(t__sim, C__i0_min);
		int__sim_max := integraal(t__sim, C__i0_max);
		int__sim := integraal(t__sim, C__i0);
		if 0 < (int__sim - int__ex)*(int__sim_max - int__ex) then
			Gamma__i_max := Gamma__i;
		else
			Gamma__i_min := Gamma__i;
		end if;
	end do;
	return Gamma__i;
end proc;

Since none of the variables are now indexed, you need to check that was what you intended - was C or C__i intended to be an Array or Vector, for example?

@Scot Gould In @acer's detailed explanation its clear that it isn't just the plot form. The operator call form is robust and is my first go-to solution for plot difficulties, especially not needing to think about where and how many uneval quotes are needed. It has been around for as long as I can remember, and the plot help page refers to it as "The commonly used operator form of the calling sequence," even if it might not be "traditional" :)  

@acer Thanks for the detailed explanation. (I realised later that just entering the GAMMA equivalent didn't give the same result, and so the real situation was more complicated.)

@SHIVAS There is an attempt to do an inverse transformation at the end of my worksheet, which shows how to do it. It does not succeed but Eq 22 was different, Perhaps if it agreed with Eq 22 it would work. You will need to find the error in the eqns or conditions (assuming there is an error) to proceed.

@SHIVAS Yes, I did not get the same answer in eqn 22. I do not know what is wrong. If you think an assumption needs to be changed you will need to say which assumption it is and what it needs to be changed to, or change it yourself.

First 20 21 22 23 24 25 26 Last Page 22 of 69