casperyc

927 Reputation

10 Badges

12 years, 23 days

MaplePrimes Activity


These are replies submitted by casperyc

Hi Carl,

Thanks! It's worked like magic!

I did not attach the code because it was a bit messy and too 'much'. And I couldnt think of an toy example. But I assume this is kind of self explanatory what my aim is.

casper

Hi Carl,

Thanks! It's worked like magic!

I did not attach the code because it was a bit messy and too 'much'. And I couldnt think of an toy example. But I assume this is kind of self explanatory what my aim is.

casper

@Alejandro Jakubi Thanks! I will see if I can identify which part of the code is slowing Maple down and takes most of the time.

@Alejandro Jakubi Thanks! I will see if I can identify which part of the code is slowing Maple down and takes most of the time.

@Carl Love Hi, higher numerical rank occured only a few times when I run batch process on the school server, in which case I did not write my program to captuer its initial (random) values or so. But ya, if I find another example or instant where it happens again, I will surely post it.

Thanks for the myRand() procedure. I will try it out!

casper

@Carl Love Hi, higher numerical rank occured only a few times when I run batch process on the school server, in which case I did not write my program to captuer its initial (random) values or so. But ya, if I find another example or instant where it happens again, I will surely post it.

Thanks for the myRand() procedure. I will try it out!

casper

I will try those ?exprofile and ?excallgraph. See what happens.

As I am a research student, I am not sure if I follow exactly what you mean by real reals, or real floats? What I did and why I did 2dp is because I want to make sure the 'p' and the 'phi' are between 0 and 1. No offense but if I put my problem simply, what I need to find is that for some complicated model, the rank follow certain pattern. The guess is the pattern 'jumps', so if possible i want to check for some 'large' values. in this case, possibly K>12.

Honestly, i dont know how to generate real reals, or real floats as well. If you see my coding, I have also tried to use linear instead of exponentials. That would give me 'some' freedom of choosing random points.

I have also had a problem, say where the symbolic rank is 4. But we get a numerical rank 5? Which I would think impossible. Numerically, it should only be smaller.

For now, as the problem of being slow, could that be evaluating

kappa:=map(indprob,CH,K,C,tvarphi,tvarp,phiMix);


As K gets larger, the number of terms increase almost like 2^K-2, so for K=12, the map does it for 4094 times, but for K=15, it has to map 32766 times. And also then the derivate matrix.

 

Really thanks for the tips and suggestion!

 

Casper

I will try those ?exprofile and ?excallgraph. See what happens.

As I am a research student, I am not sure if I follow exactly what you mean by real reals, or real floats? What I did and why I did 2dp is because I want to make sure the 'p' and the 'phi' are between 0 and 1. No offense but if I put my problem simply, what I need to find is that for some complicated model, the rank follow certain pattern. The guess is the pattern 'jumps', so if possible i want to check for some 'large' values. in this case, possibly K>12.

Honestly, i dont know how to generate real reals, or real floats as well. If you see my coding, I have also tried to use linear instead of exponentials. That would give me 'some' freedom of choosing random points.

I have also had a problem, say where the symbolic rank is 4. But we get a numerical rank 5? Which I would think impossible. Numerically, it should only be smaller.

For now, as the problem of being slow, could that be evaluating

kappa:=map(indprob,CH,K,C,tvarphi,tvarp,phiMix);


As K gets larger, the number of terms increase almost like 2^K-2, so for K=12, the map does it for 4094 times, but for K=15, it has to map 32766 times. And also then the derivate matrix.

 

Really thanks for the tips and suggestion!

 

Casper

Bug?

Why is it working in the "working" example then? Cost me hours already!

Thanks for the get around!

Bug?

Why is it working in the "working" example then? Cost me hours already!

Thanks for the get around!

@acer 

 

restart:

with(LinearAlgebra):

interface(rtablesize=infinity);

intM:=proc(K::integer, C::integer)

description "Creat interaction terms (matrix form).";

local
    M:=Matrix(K,C,symbol=mix),
    tolt,cols,rows;

# Set constraints

    M(..,1):=0;
    M(1..2,..):=0;
    
    tolt:=convert(M,`+`);
    cols:=MTM:-sum(M,1);
    rows:=MTM:-sum(M,2);
    
    M(2,..):=-cols;
    M(..,1):=-rows;
    M(2,1):=tolt;
    
    M;
    
end proc;


indprob:=proc(x,
    K::integer,
    C::integer,
    tvarphi::truefalse:=false,
    tvarp::truefalse:=false,
    phiMix::truefalse:=false
    )

description "Compute the Probability for given capture history $x$.";

local
    phi,p, # Vectors of phi and p
    int0,  # Matrix of interaction terms
    ind,   # Find indices of where the '1's are
    fcap,  # First capture
    li,    # Last capture
    ans;   # Desired probability

# global phi,p; # Vectors of phi and p;

# take care of phi (nu)
# SINGLE CLASS
    
    if tvarphi then
        phi:=Array(1..K,[seq(1/(1+exp(-nu[j])),j=1..K)]);
        # phi:=Array(1..K,[seq(nu[j],j=1..K)]);
    else
        phi:=Array(1..K,[seq(1/(1+exp(-nu)),j=1..K)]);
        # phi:=Array(1..K,[seq(nu,j=1..K)]);
    end if;
    phi(K):=0;
    
# take care of p

    if tvarp then
        if phiMix then
            int0:=intM(K,C);
            p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+tau[j]+eta[c]+int0[j,c]))));
            # p:=Array(1..K,1..C,(j,c)->mu+tau[j]+eta[c]+int0[j,c]);
        else
            if C=1 then
                p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+tau[j]))));
                # p:=Array(1..K,1..C,(j,c)->mu+tau[j]);
            else
                p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+tau[j]+eta[c]))));
                # p:=Array(1..K,1..C,(j,c)->mu+tau[j]+eta[c]);
            end if;
        end if;
    else
            p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+eta[c]))));
            # p:=Array(1..K,1..C,(j,c)->mu+eta[c]);
    end if;
    
    p(1,..):=0;
    p:=subs({tau[2]=0,eta[1]=0},p);

    ind:=ArrayTools[SearchArray](Array(x));
    fcap,li:=ind[1],ind[-1];
    
    ans:=
        add(
            w[c]*
            add(
                (1-phi[d])*
                mul(phi[j],j=fcap..(d-1))
                *
                mul(p[j,c]^(x[j])*(1-p[j,c])^(1-x[j]),j=(fcap+1)..d),
                d=li..nops(x)
            ),
        c=1..C);
    
    ans;
    
end proc;


myP:=proc(
    K::integer,
    C::integer,
    tvarphi::truefalse:=false,
    tvarp::truefalse:=false,
    phiMix::truefalse:=false,
    n::integer:=30)

description "Test for Parameter Redundancy (for mixtures of $p$).";    

local
    CH,      # All Capture Histories from $K$ sample
    kappa,   # Exhaustive Summary
    pars,    # ALL Parameters in the model
    DD1,     # Derivative Matrix (Jacobian) of Kappa
    
    ## Numerical ##
    
    ms,      # Dummy var, loop
    myrand,  # Random number generator
    numpars, # Initial parameters
    D1rand,  # Matrix with numbers
    rNtmp,   # Rank from D1rand
    rS,      # Symbolic Rank
    rN       # Numerical Rank
    ;
    
# Compute Capture Histories

    CH:=combinat[permute]([seq(s$K,s=[0,1])],K);
    CH:=CH[2..-1];

# Get Exhaustive Summary

    kappa:=map(indprob,CH,K,C,tvarphi,tvarp,phiMix);
    kappa:=subs({w[C]=1-add(w[i],i=1..(C-1))},kappa);
    kappa:=subs({seq(w[i]=1/(1+exp(-wnew[i])),i=1..(C-1))},kappa);
    kappa:=kappa[2..-1];
    
# Get Jacobian
    pars:=indets(kappa,name);
    DD1:=VectorCalculus[Jacobian](Vector(kappa),convert(pars,list));

    #rS:=Rank(DD1);

    rN:=0;

    myrand:=evalf(rand(-300..300)/100);

    for ms to n do

        numpars:=seq(pars[i] = myrand(), i = 1 .. nops(pars));

        D1rand:=eval(DD1, {numpars});

        #D1rand:=evalf[3](D1rand);
        D1rand:=evalf(D1rand);

         rNtmp:= LinearAlgebra[Rank](D1rand);

         if (rNtmp>rN) then
             rN:=rNtmp;
         end if;

    end do;
    
    return rN;
    
end proc;


################################################
################################################
myCK:=proc(kmin::integer,kmax::integer,C::integer,n)

##################################################
    print(seq(myP(K,1,false,false,false,n),K=kmin..kmax));
    print(seq(myP(K,1,false,true,false,n),K=kmin..kmax));
    print(seq(myP(K,C,false,false,false,n),K=kmin..kmax));
    print(seq(myP(K,C,false,true,false,n),K=kmin..kmax));
    print(seq(myP(K,C,false,true,true,n),K=kmin..kmax));
##################################################
printf("\nNext block...\n");
    print(seq(myP(K,1,true,false,false,n),K=kmin..kmax));
    print(seq(myP(K,1,true,true,false,n),K=kmin..kmax));
    print(seq(myP(K,C,true,false,false,n),K=kmin..kmax));
    print(seq(myP(K,C,true,true,false,n),K=kmin..kmax));
    print(seq(myP(K,C,true,true,true,n),K=kmin..kmax));
##################################################
printf("Done...\n");
end proc;


And some results that I currectly have:









I wonder if it is possible to speed it up ? As I can't really get any result for ''larger model" as in

myCK(3,BIG NUMBER>12,C,10);

I "think" increasing C is "OKish", but the real problem is increasing K>10 (2nd input). It's really slow and most of the time, I just dont have an answer for days...

 

Many thanks!

 

casper

What's 'real time' and CPU time exactly? Which one says "efficient" (or "quicker")? I tried to look them up but I still dont quite get it. Clearly from your example, one tells Map is better and the other tells map is better.

What's 'real time' and CPU time exactly? Which one says "efficient" (or "quicker")? I tried to look them up but I still dont quite get it. Clearly from your example, one tells Map is better and the other tells map is better.

And there was a mistake in my post. I did try

map2(myfunc, phi, p1, {CH||(1..5)});

and got a error with

"dimensions must be specified with this type of initializer".

 

Yes, i will look into the Threads library!

Never tried it.

 

Is it not by default for Maple to 'adapt' to computer and use multi thread then?

Thanks,

Casper

And there was a mistake in my post. I did try

map2(myfunc, phi, p1, {CH||(1..5)});

and got a error with

"dimensions must be specified with this type of initializer".

 

Yes, i will look into the Threads library!

Never tried it.

 

Is it not by default for Maple to 'adapt' to computer and use multi thread then?

Thanks,

Casper

First 7 8 9 10 11 12 13 Last Page 9 of 22