Joe Riel

9660 Reputation

23 Badges

20 years, 3 days

MaplePrimes Activity


These are answers submitted by Joe Riel

If I understand correctly, there isn't a way to do what you want. That is, given an existing module there is no straightforward way to add an export to that module. Even modifying an export of a module isn't generally possible (there are some hackish ways, but they have limits). This is a drawback of the modular approach compared to the old table-based package approach used in Maple R5—while the modular approach is more robust and better from the programmar's (creator's) viewpoint, it is less malleable for the user. It is, of course, possible to add a new module to an existing archive, though usually it is better to just create a new archive and insert it in the libname path. Use LibraryTools[Save]. To make a module ``withable,'' give it the package option. For example:
mymodule := module()
export A,B:
option package;
  A := proc() ... end proc:
  B := proc() ... end proc:
end module:
LibraryTools:-Save(mymodule,"mylib.mla"):
The ability to rebind operators, such as +, was introduced in Maple 6; however, the overload option (and procedure) was introduced in Maple 9.5.
See the help page for NextAfter.
> Digits := 20:             
> epsilon := NextAfter(0,1);
                                                                   -20
                                                  epsilon := 0.1 10

> 1+epsilon;                
                                                 1.0000000000000000000
That isn't quite what you want, but...
I'm not sure I understand, but the following (rather ugly) code does what I think you want
restart; with(stats); Digits := 14;

N := 4:
L := 0.002:
beta:=0.85:
cost := Vector([200,100,700,1000]):
lambda := Vector([0.2,0.6,0.3,0.7]):
Lambda := add(lambda[i], i=1..N):

DFR := s -> ( add(lambda[i]*add((L*lambda[i])^k*exp(-L*lambda[i])/k!
                                , k = 0 .. s[i]-1)
                  , i = 1 .. N)
              / Lambda );

IN := Matrix(N,N,shape=identity):
S := Vector[row](4):
maxval := DFR(S);

while maxval < beta do
    prevval := maxval;
    maxrelval := 0;
    for k to N do
        Sk := S + IN[k,1..N];
        val := DFR(Sk);
        relval := (val - prevval)/cost[k];
        if relval > maxrelval then
            maxrelval := relval;
            maxval := val;
            Snxt := Sk;
        end if;
    end do;
    if maxrelval = 0 then
        print("cannot increase");
        break;
    end if;
    S := Snxt;
end do:

S, maxval;
I think the following does what you want.
restart; with(stats); Digits := 14;

N := 4:
L := 0.002:
beta:=0.85:
lambda := Vector([0.2,0.6,0.3,0.7]):
Lambda := add(lambda[i], i=1..N):

DFR := s -> ( add(lambda[i]*add((L*lambda[i])^k*exp(-L*lambda[i])/k!
                                , k = 0 .. s[i]-1)
                  , i = 1 .. N)
              / Lambda );

IN := Matrix(N,N,shape=identity):
S := Vector[row](4):
maxval := DFR(S);

while maxval < beta do
    found := false; # flag to avoid infinite loop if no Sk is better than S.
    for k to N do
        Sk := S + IN[k,1..N];
        val := DFR(Sk);
        if val > maxval then
            maxval := val;
            Snxt := Sk;
            found := true;
        end if;
    end do;
    if not found then break; end if; # may want to add warning
    S := Snxt;
end do;

maxval;
I don't understand what you are attempting, however, there seem to be two clear problems. First, you probably want to use add instead of sum in the computation of DFR (you used one of each). Second, your do loop is confusing. It's index is i, but you also use i as the index in the adds. These are different i's and I'm not sure which you intend. Change one of them to a j to make it clear what you are doing.
A good method for creating types is to use the TypeTools module, in particular, its AddType export. You'd do something like:
Automatism := module()
export zeros;
local ModuleLoad,ModuleUnload;

  ModuleLoad := proc()
    TypeTools:-AddType( 'TransferFunction', handler );
  end proc;

  ModuleUnload := proc()
    TypeTools:-RemoveType( 'TransferFunction' );
  end proc;

  zeros := proc( f:: TransferFunction) ... end proc;

  ModuleLoad();
end module:
You have to write the 'handler'; see the AddType help pages for details (it wasn't clear to me what your data structure is, otherwise I'd have suggested something). ,
You'd use seq something like this:
display([ seq( plot([X1[k],Y1[k]],[X2[k],Y2[k]]), k=1..20) ], 'insequence');
Note that your typing of v is not valid, try executing the procedure without compilation and kernelopts(assertlevel=2). Also, ?Compiler,Compile states:
Array Support
- Arrays can only appear as parameters of a procedure that is to be compiled.
While not ideal, you could make v a parameter:
kernelopts(assertlevel=2);
                                       0

v := Vector([1$4],datatype=float[8]):
foo := proc(x::float, v::Vector(4,datatype=float[8]))
    return x*v[2];
end proc:

foo(1.0,v);
                                      1.0


fooc := Compiler:-Compile(foo);
fooc := proc()
option call_external, define_external(_ma78ea510b1b4af6a145704a5af10ab6b,
MAPLE, LIB = "/var/tmp/joe-30090/_ma78ea510b1b4af6a145704a5af10ab6b.so");
    call_external(0, 1073846032, true, args)
end proc


fooc(1.0,v);
                                      1.
You can also use a regular conditional: SKIPLINES := true:
if SKIPLINES <> true then
...
end if;
to skip lines of code.
One way is to symbolically evaluate the function to generate an expression, then use unapply to create a new function:
g := unapply(f(A,'b'),'b');
                              g := b -> 3.5 + b
Note the forward quotes around b to delay its evaluation. This method won't always work, since it isn't always possible to symbolically evaluate a function. A more general technique is the following:
g := subs(_A=A, b -> f(_A,b));
                                   g := b -> f(3.5, b)
Another way to do that is the following
g := (a -> b -> f(a,b))(A);
                                   g := b -> f(3.5, b)
It is briefly mentioned in the the help page for rtable. Its first argument is a range (0..1), its second argument specifies the percentage of nonzero values (1.0, you want all to be nonzero). The elements of the rtable are randomly and uniformly initialized between 0 and 1. Using frandom is much faster than other methods (see generating an array of random floats). Here is an example of creating a Vector with 10 random elements:
rnd := rtable(1..10, frandom(0..1,1.0), subtype=Vector[column], datatype=float[8]):
                                                  [ 0.530484209985166988 ]
                                                  [ 0.952299352583890979 ]
                                                  [0.0354787515261802966 ]
                                                  [ 0.240688559954775894 ]
                                                  [ 0.767101828106371397 ]
                                         rnd :=   [ 0.724353985897650987 ]
                                                  [ 0.388089818893285866 ]
                                                  [ 0.219477841874477498 ]
                                                  [0.00149287817728582796]
                                                  [ 0.436224435569102509 ]
You can use the frandom initializer to rtable to quickly create a Vector of floats. In the following I'll generate a Matrix (rather than a Vector) so that I can use the PlotHistogram procedure from ImageTools to generate the histogram.
rnd := rtable(1..10\000, 1..1, frandom(0..1,1.0), subtype=Matrix, datatype=float[8]):
plt := ImageTools:-PlotHistogram(rnd,10):
display[plots](plt, style=point);
Just do ans mod d*e. Actually, because the inline mod operator can be reassigned to be either modp or mods (see the help page for mod), it is frequently better to use the modp or mods function directly. Thus, change the last line in the procedure to
modp({...}, d*e)
, where the ellipsis indicates the original statement (depending on which version you used.
Good point. Note that the problem exists even if you stay in document mode (or 2D input). That is, cutting and pasting the output to an input region causes a conversion from the literal subscript format to the indexed format. Rather than doing evalb(%=t[0]), you can evaluate lprint(%), that prints the 1D representation of the symbol. In this case,
lprint(t0);
`#msub(mi("t"),mn("0"))`
After pasting the output into an lprint we get
lprint(t0);
t[0]
I rarely cut and paste from the output into an input region. Better to learn to use Maple operators to pick off the element of interest.
First 104 105 106 107 108 109 110 Last Page 106 of 114