Carl Love

Carl Love

28015 Reputation

25 Badges

12 years, 297 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@lcz And here is another use of RelationGraph: a procedure for GraphPower with 3 methods:

  1. an algorithm that I just invented called unions;
  2. a method based on the AllPairsDistance matrix and RelationGraph;
  3. Maple's own GraphTheory:-GraphPower.

My unions algorithm in several times faster than the others for small values of k.

RelationGraph:= (f, V::list({integer, symbol, string}))->
local n:= {$nops(V)}, i, F:= i-> j-> f(V[i], V[j], _rest);
    GraphTheory:-Graph(V, [seq](select(F(i), n, _rest), i= n))
:

GraphPower:= proc(
    G::Graph, 
    k::posint,
    {method::identical("default", "distance", "unions"):= ""},
    $
)
local V:= GraphTheory:-Vertices(G), n:= [$nops(V)], N, L;
    if method = "unions" then
        L:= rtable(op~((N:= op(4,G))));
        to k-1 do N:= N union~ (n-> {seq}(L[[n[]]]))~(N) od;
        GraphTheory:-Graph(V, N minus~ `{}`~(n))
    elif method = "distance" then
        N:= GraphTheory:-AllPairsDistance(G);
        L:= table(V=~ n);
        RelationGraph((u,v)-> u<>v and N[L[u],L[v]] <= k, V)
    else
        GraphTheory:-GraphPower(G,k)
    fi    
end proc
:

 

And here is a usage test comparing all methods:

GT:= GraphTheory:  RG:= GT:-RandomGraphs:
G:= RG:-RandomGraph(2^10, 2^11):
gc(): CodeTools:-Usage(GraphPower(G, 5, method= "default")):
memory used=71.82MiB, alloc change=6.01MiB, 
cpu time=1.16s, real time=1.15s, gc time=0ns

E1:= GT:-Edges(%):  nops(%);
                             315735

gc(): CodeTools:-Usage(GraphPower(G, 5, method= "unions")):
memory used=60.46MiB, alloc change=0 bytes,
cpu time=281.00ms, real time=280.00ms, gc time=0ns

E2:= GT:-Edges(%):  evalb(E1=E2);
                              true

gc(): CodeTools:-Usage(GraphPower(G, 5, method= "distance")):
memory used=240.44MiB, alloc change=6.01MiB, 
cpu time=3.69s, real time=2.99s, gc time=812.50ms

E3:= GT:-Edges(%):  evalb(E1=E3);
                              true

 

@lcz You wrote:

  • I am probably more interested in the panning implementation of RelationGraph in maple. Maybe it's not always the most effective. But it can reflect the degree of abstraction of certain graph operations. And not just some graph operations in isolation.

I agree on all points. A procedure for RelationGraph is quite easy; here it is:

RelationGraph:= (f, V::list({integer, symbol, string}))->
local n:= {$nops(V)}, i, F:= i-> j-> f(V[i],V[j]);
    GraphTheory:-Graph(V, [seq](select(F(i), n, _rest), i= n))
:


And here is a rewrite of LexProd which offers a choice of two methods: one using RelationGraph and the other the same as the algorithm in my previous Reply:

LexProd:= proc(
    G::seq(Graph),
    {
        method::identical("default", "relation"):= "default",
        vertex_format::{string, symbol}:= "(%a,%a)"
    },
    $
)
local
    LexProd2:= proc(G::Graph, H::Graph, $)
    local 
        (Vg,Vh):= op(GT:-Vertices~([G,H])),
        (ng,nh):= op(`..`~(1, nops~([Vg,Vh]))),
        (Ng,Nh):= op(op~(4, [G,H])),
        i, j, J, k, K,
        P:= [seq](seq([i,j], j= nh), i= ng),
        Vgh:= (curry(nprintf, vertex_format)@((i,j)-> (Vg[i],Vh[j]))@op)~(P)
    ;
        if method="relation" then
            K:= subs(_K= op~(table(Vgh=~ P)), v-> _K[v]);
            RelationGraph(
                proc(a, b, $)
                local (u,v):= K(a), (x,y):= K(b); 
                    x in Ng[u] or u=x and y in Nh[v]
                end proc,
                Vgh
            )
        else
            k:= 0;  K:= table((p-> op(p)= ++k)~(P));
            GraphTheory:-Graph(
                Vgh,
                [seq](
                    {seq(seq(K[k,J], k= Ng[j[1]]), J= nh), seq(K[j[1],k], k= Nh[j[2]])},
                    j= P             
                )  
            )
        fi
    end proc,
    n:= nops([G]);
;
    if n=0 then error "no graphs given" elif n=1 then G[1] else foldl(LexProd2, G) fi
end proc
:


And here's an example verifying that the two methods give identical results:

GT:= GraphTheory:
(G,H):= GT:-Graph~([`$`]~([3,4]), [{{1,2}}, {{1,3},{2,4}}])[];
 G, H := Graph 38: an undirected unweighted graph with 3 vertices and 1 edge(s),
     Graph 39: an undirected unweighted graph with 4 vertices and 2 edge(s)

GH1:= LexProd(G,H);
    GH1 := Graph 40: an undirected unweighted graph with 12 vertices and 22 edge(s)

GH2:= LexProd(G, H, method= "relation");
    GH2 := Graph 42: an undirected unweighted graph with 12 vertices and 22 edge(s)

evalb(GT:-Edges(GH1) = GT:-Edges(GH2));
                              true

 

@Earl You wrote:

  • I've spent some time studying your latest replay. It's difficult for me to understand so I'll return to it later.

I'd be happy to answer any questions about it, especially if it's not a large number of questions asked all at once. (However, a large number of questions spread out over several sessions is not a problem.)

@Joe Riel @AmirHosein Sadeghimanesh

It should be pointed out, for the OP's benefit, that

  1. Joe's workaround works because MutableSet has a procedure named ModuleIterator;
  2. the original command doesn't work because MutableSet doesn't have a procedure named ModuleType;
  3. the workaround in my Answer below works because MutableSet does have a procedure named convert which accepts the keyword set as its 2nd argument.

In summary, what works is determined entirely by an object's defining module having specifically named procedures (or "methods" if you want to call procedures that).

@planetmknzm

Yesterday I saw that you asked (essentially) how to compose an animation consisting of some parts that move and some that don't. That Question got deleted, I don't know if by you or some Moderator. Anyway, I'll Answer that Question here.

Here is a program that plots a curve on the unit sphere and animates a plane tangent to the sphere whose point of tangency moves along the curve. As you can see, this can be done much more simply than your attempts. My computer uses 0.17 seconds to produce this 100-frame animation.

restart:
Sph:= (theta,phi)-> [(cos,sin)(theta)*~sin(phi), cos(phi)]:
C:= Sph@op@[theta, phi]:

#Example curve:
theta:= t-> t+Pi:  phi:= t-> t^2/Pi:

#Parts of the animation that don't move:
Static:= plots:-display(
    plottools:-sphere(), 
    plots:-spacecurve(C(t), t= -Pi..Pi, thickness= 4, color= blue)
):
#Plot a plane tangent to sphere at a given point P: 
Plane:= proc(P)
local v:= <P>^%T, B:= LinearAlgebra:-NullSpace(v)^~%T;
    plots:-polygonplot3d(<v+B[1], v+B[2], v-B[1], v-B[2]>)
end proc
: 
plots:-animate(Plane, [C(t)], t= -Pi..Pi, frames= 100, background= Static);

@Joe Riel The OP wants to check the type of the elements of the MutableSet. So, continuing with your example,

M:= MutableSet(a, b, c);
type(M, 'MutableSet'(symbol));

Error, module does not have a ModuleType member to accept structured type arguments

Please give a simple example of esp for which it doesn't work. As far as I can tell, what you decribed works for me:

esp:= 2*x + 3*x^2 + 4*sqrt(x):
coeff(esp, sqrt(x));

                               
4

@Jno What is your Maple version (e.g., I'm using Maple 2022)?

@Jno 

I've updated the worksheet in my previous Reply. I think that it makes the math easier to understand, although the Maple coding is a bit trickier. In particular, I came up with a single simple unified formula for combining sums of squares and their degrees of freedom, which simplifies the presentation of SS_ASS_AP, and SSE. Please download the new worksheet.

The "error" (better termed the "unmodeled variance") is a factor of all the F-ratios, so it would have been impossible to proceed without it.

Maple's Vector and Matrix indices always start at 1. Array indices can start at any integer. I think that you're misinterpretting the DataFrame's column of row labels and row of column labels as being a 0th column and 0th row; they're not---they're just labels shown as a visual aid. They're akin to the tickmarks on a plot: a visual aid that's not part of the function being plotted.

Negative Vector and Matrix indices represent counting rows or columns backwards, i.e., from the bottom or right ends. Zero is never a valid Vector or Matrix index.

@AdeWaele What exactly do you find "not easily possible" about acer's Answer dgdx:=  D[1](g)??? It's the other two Answers that are flawed!

@Jno 

Re HFloat(...): That stands for "hardware float", which is the same as "double precision" in several other languages. It's correct, but unneeded information. I'm not sure why you're seeing those---perhaps  a display setting. Nonetheless:

  1. There's no reason that you need to see them.
  2. They can be ignored.
  3. The following program update removes them anyway.

Re accessing/extracting the values: I've changed the output format to a DataFrame (very similar to a Matrix) that presents a traditional ANOVA table. At the bottom of the following worksheet, I show how to extract/access the values.

I extended the procedure to do the complete two-way ANOVA with F-ratios and their p-values.

restart
:

ANOVA_2way:= proc(X::And(Array, 3 &under ArrayNumDims))
description `Two-factor analysis of variance (aka, two-way ANOVA)`;
option `Author: Carl Love <carl.j.love@gmail.com> 2022-Jul-1`;
uses St= Statistics;
local
    AAD:= (A::rtable, d::posint)-> #replacement for ArrayTools:-AddAlongDimension
    local k;
        add(A[(..)$d-1, k], k= [ArrayDims](A)[d]),

    N:= numelems(X), #total number of observations
    Xij:= AAD(X, 3), #sums of trials  

    #sum-of-squares procedures:
    SS:= (A::rtable)-> local n:= numelems(A); (add(A^~2)*n/N, n),
    all:= SS(<add(Xij)>)[1], #ie, Sum(X)^2/N,
    SSrec:= (SS::realcons, n::And(posint, Not(1)))->
        Record("df"= n-1, ("SS","MS")=~ (SS-all)/~(1, n-1)),
    `&-`:= (R::record, S::record)->
        SSrec(R:-SS - S:-SS + all, R:-df - S:-df + 1),

    #sums of squares due to...
    SS_P:= SSrec(SS(AAD(Xij, 2))), #...partition ("P factor")
    SS_A:= SSrec(SS(AAD(Xij, 1))), #...treatment ("A factor")
    SS_AP:= SSrec(SS(Xij)) &- SS_P &- SS_A, #...interaction (P x A)
    TSS:= SSrec(SS(X)), #...total, ie, Sum(X^~2) - Sum(X)^2/N
    SSE:= TSS &- SS_P &- SS_A &- SS_AP, #...random or unaccounted factors {"error")

    #procedure for each row of output (with F-ratio and its p-value):
    F:= (SS::record)->
    local F:= SS:-MS/SSE:-MS;
        <
            SS:-SS | SS:-df | SS:-MS | F |
            1 - St:-CDF('FRatio'(SS:-df, SSE:-df), F, 'numeric')
        >,

    Out:= <F~([SS_P, SS_A, SS_AP, SSE, TSS])[]> #ANOVA table
;
    #Put blanks for N/A entries at bottom right:
    Out[-2, -2..]:= ``;  Out[-1, -3..]:= ``;
    DataFrame(
        #Change HFloats and round to 6 sig. digits:
        subsindets(Out, float, evalf[6]),         
        'rows'= ["P", "A", "PxA", "Error", "Total"],
        'columns'= ["SS", "df", "MSS", "F", "p-val"]
    )
end proc
:

#Example usage:

#Appraiser (A) - Part (P) data with 3 measurements for each (A,P) pair:
AP:= Array(
    [
         [[29, 41, 64],    [8,   25,  7],    [4, -11, -15]],
        -[[56, 68, 58],    [47, 122, 68],    [138, 113, 96]],
         [[134, 117, 127], [119, 94, 134],   [88, 109, 67]],
         [[47, 50, 64],    [1, 103, 20],     [14, 20, 11]],
        -[[80, 92, 84],    [56, 120, 128],   [146, 107, 145]],
         [[2, -11, -21],   [-20, 22, 6],    -[29, 67, 49]],
         [[59, 75, 66],    [47, 55, 83],     [2, 1, 21]],
         [-[31, 20, 17],   [-63, 8,- 34],   -[46, 56, 49]],
         [[226, 199, 201], [180, 212, 219],  [177, 145, 187]],
        -[[136, 125, 131], [168, 162, 150],  [149, 177, 216]]
    ]/100,
    datatype= hfloat
):

R:= ANOVA_2way(AP);

"DataFrame([[[88.3619,9,9.81799,213.517,0.],[3.16726,2,1.58363,34.4401,1.09385 10^(-10)],[0.358982,18,0.0199435,0.433721,0.974106],[2.75893,60,0.0459822,,],[94.6471,89,,,]]],rows=["P","A","PxA","Error","Total"],columns=["SS","df","MSS","F","p-val"])"

#Access/extract the values like this:
R["A", "SS"];
R["PxA", "p-val"];

3.16726

.974106

#...or matrix style, like this:
R[2,1];
R[3,5];

3.16726

.974106

 

Download ANOVA.mw

@C_R Hmm, that link worked for me earlier, but it doesn't anymore. It was a 3-4-page paper with 2 examples of solving a 7th-order bivariate PDE by the Variational Iteration Method. The examples were done in Maple, but the commands weren't explicitly shown.

@Jno You don't have any Array named TestA. If TestA existed, then if it caused an error, the error message would list the contents of TestA rather than the name "TestA".

The command restart erases all variables from the "active" memory. Thus, if you do something like the following, you'll get that error:

TestA:= Array([1]);
restart;
ArrayDims(TestA);

A correct order is

restart;
TestA:= Array([1]);
ArrayDims(TestA);

@mehran rajabi Please post an executed worksheet that shows the blank space where the plot should be.

@Jno My guess is that you declared nm1 as an array rather than as an Array. The lowercase array has been obsolete for about 20 years. It only still exists so that old code will still run. Likewise, vector and matrix are obsolete.

Arrays, Vectors, and Matrixes are collectively called rtables. My line of code that you asked about specifies that the procedure ANOVA_3D has one parameter, X, whose corresponding passed argument must be both an rtable And have exactly 3 dimensions.

First 85 86 87 88 89 90 91 Last Page 87 of 708