apost3a.html

mafapo3c.mws

Depuração e geração de programas C

Routo Terada

Depuração (Debugging)

> restart;

> gera:=rand(1..5);

> gera();

> printlevel:=25; # definir nível de info exibida
# para auxiliar debugging

> gera();gera();

{--> enter gera, args =

<-- exit gera (now at top level) = 1}

{--> enter gera, args =

<-- exit gera (now at top level) = 3}

> printlevel:=1; # volta para o nível 1 (default)

> gera();

> funfat:=proc(n) # função fatorial de n = n*(n-1)!
if n<2 then 1;
else n*funfat(n-1);
fi;
end;

> printlevel:=25;

> funfat(4);

{--> enter funfat, args = 4

{--> enter funfat, args = 3

{--> enter funfat, args = 2

{--> enter funfat, args = 1

<-- exit funfat (now in funfat) = 1}

<-- exit funfat (now in funfat) = 2}

<-- exit funfat (now in funfat) = 6}

<-- exit funfat (now at top level) = 24}

>

> printlevel:=1;

Geração de código C com codegen

> with(codegen):

Warning, the protected name MathML has been redefined and unprotected

>

> codegen[C](funfat);

#include <math.h>

double funfat(n)

double n;

    if( n < 2.0 )

      return(1.0);

    else

      return(n*funfat(n-1.0));

> f := 1-x/2+3*x^2-x^3+x^4;

> codegen[C](f);

      t0 = 1.0-x/2.0+3.0*x*x-x*x*x+x*x*x*x;

> codegen[C](f,optimized);

      t2 = x*x;

      t5 = t2*t2;

      t6 = 1.0-x/2.0+3.0*t2-t2*x+t5;

> f := Pi*ln(x^2)-sqrt(2)*ln(x^2)^2;

> codegen[C](f,optimized);

      t1 = x*x;

      t2 = log(t1);

      t4 = sqrt(2.0);

      t5 = t2*t2;

      t7 = 0.3141592653589793E1*t2-t4*t5;

> cs := [s = 1+x, t = ln(s)*exp(-x), r = exp(-x)+x*t ];

> codegen[C]( cs, optimized );

      s = 1.0+x;

      t1 = log(s);

      t2 = exp(-x);

      t = t1*t2;

      r = t2+x*t;

> v := array([exp(-x)*x,exp(-x)*x^2]);
codegen[C](v,optimized);

      t1 = exp(-x);

      t3 = x*x;

      v[0] = t1*x;

      v[1] = t1*t3;

> A := array(1..2,1..2,symmetric):
A[1,1] := log(x): A[1,2] := 1-log(x):
print(A);

> codegen[C](A,optimized); # Note um elemento indefinido

      t1 = log(x);

      t2 = 1.0-t1;

      A[0][0] = t1;

      A[0][1] = t2;

      A[1][0] = t2;

      A[1][1] = (0.0/0.0);

>

Geração de código otimizado

Calcular 1+y+y^2/2+...+y^n/n! como uma soma simbólica, onde y = (1-x)^2.

O tradutor C não conhece o comando sum.

> with(codegen,prep2trans):
g := proc(x,n)
local i,t;
t := (1-x)^2;
sum(t^i/i!,i=0..n);
end:
codegen[C](g);

#include <math.h>

double g(x,n)

double x;

double n;

  double i;

  double t;

    t = pow(1.0-x,2.0);

codegen/C/expression:   Unknown function:   sum   will be left as is.

codegen/C/expression:   Unable to translate   i = 0 .. n   into C

codegen/C/expression:   Replacing value with three parameters   i, 0, n

    return(sum(pow(t,1.0*i)/exp(lgamma(1.0+i)),i,0.0,n));

Abaixo, com o codegen, um for é gerado correspondendo ao comando sum

> g := prep2trans(g,language=C);

> codegen[C](g);

double g(x,n)

double x;

double n;

  double i;

  int i2;

  double s2;

  double t;

  double t1;

  double t2;

    t1 = 1.0-x;

    t = t1*t1;

    if( n < 0.0 )

      s2 = 0.0;

    else

        t2 = 1.0;

        s2 = 1.0;

        for(i2 = 1;i2 <= n;i2++)

          t2 = t/i2*t2;

          s2 += t2;

    return(s2);

A seguir um proc Maple que aceita um array unidimensional como entrada.

Nesta versão, o tradutor não conhece os limites do array.

Na versão seguinte informaremos os limites com declarations .

> h := proc(n,a::array,x) local i,s;
s := a[n];
for i from n-1 by -1 to 0 do s := a[i]+x*s; od;
s
end:

codegen[C](h);

codegen/C/pretreatment:   WARNING: Unable to determine array bounds in   a[n]

codegen/C/pretreatment:   WARNING: Unable to determine array bounds in   a[i]

double h(n,a,x)

int n;

double *a;

double x;

  int i;

  double s;

    s = a[n-1];

    for(i = n-1;0.0 <= i;i--)

      s = a[i-1]+x*s;

    return(s);

> codegen[C](h,declarations='a::array(1..n)');

double h(n,a,x)

int n;

double *a;

double x;

  int i;

  double s;

    s = a[n-1];

    for(i = n-1;0.0 <= i;i--)

      s = a[i-1]+x*s;

    return(s);

A seguir um exemplo de proc que cria e retorna uma matriz de tamanho constante

> g := proc(x,y,z)
local dt,grd,t;
grd := array(1 .. 3);
dt := array(1 .. 3);
dt[3] := 2*z;
t := z^2;
grd[1] := cos(x)*z-sin(x)*t;
grd[2] := 0;
grd[3] := sin(x)+cos(x)*dt[3]-1/t^2*dt[3];
grd
end:
codegen[C](g);

#include <math.h>

void g(x,y,z,grd)

double x;

double y;

double z;

double grd[3];

  double dt[3];

  double t;

    dt[2] = 2.0*z;

    t = z*z;

    grd[0] = cos(x)*z-sin(x)*t;

    grd[1] = 0.0;

    grd[2] = sin(x)+cos(x)*dt[2]-1/(t*t)*dt[2];

    return;

>

The code generation package codegen

Calling Sequence:

codegen[function](args)

Description:

The code generation package is a collection of tools for creating, manipulating, and translating Maple procedures into other languages. This includes tools for automatic differentiation of Maple procedures, code optimization, translation into C and Fortran, and an operation count of a Maple procedure.
The functions available in the package are:

C fortran optimize cost
GRADIENT HESSIAN JACOBIAN horner
makeglobal makeparam makeproc makevoid
dontreturn packargs packlocals swapargs
renamevar declare intrep2maple maple2intrep
prep2trans split joinprocs

For more information on a particular function see codegen[function] . See also codegen[example] for a worked application example.

To use a codegen function, either define that function alone using the command with(codegen, function) , or define all codegen functions using the command with(codegen) . Alternatively, invoke the function using the long form codegen[function] . This long form notation is necessary whenever there is a conflict between a package function name and another function used in the same session.

Estatística

Routo Terada

Método dos mínimos quadrados: ajuste por curva

> with(stats[fit]); # carrega o package de ajuste

> ValoresX:=[1,2,3,4];

> ValoresY:=[1,5,24.5,52.3];

Agora aplicamos o método para calcular os coeficientes de

> leastsquare[ [x, y], y=a*x^2+b*x+c, {a, b, c} ]([ValoresX, ValoresY]);

> plot([5.95*x^2-12.41*x+7.1,[[1,1],[2,5],[3,24.5],[4,52.3]]],x=1..4,linestyle=[2,3],color=[red,blue],thickness=3);

Exemplo 2 - a seguir outro exemplo de fit(), com inteiros

> leastsquare[[x,y,z]]([[1,2,3,5],[2,4,6,8],[3,5,7,10]]);

Exemplo 3 - a seguir outro exemplo de fit(), também com inteiros, mas com cálculo de resíduos

> valX:=[1,2,3,4];

> valY:=[0,14,6,24];

> eq_fit:= leastsquare[[x,y], y=a*x^2+b*x+c]([valX, valY]);

Vamos traçar as duas curvas a seguir

> tab:=[]; j:=1;
while j<=nops(valX) do
tab:= [ op(tab), [valX[j], valY[j]] ];
j:=j+1;
od;

> plot([tab, x^2+7/5*x],x=0..5, color=[red,blue], style=[line,line],thickness=2);

Como transformar eq_fit em procedure, para calcular o resíduo depois

> eq_func:=unapply(rhs(eq_fit),x);

> with(stats);

Warning, these names have been redefined: anova, describe, fit, importdata, random, statevalf, statplots, transform

Depois, calculamos os valores previstos:

> valYprevisto:=transform[apply[eq_func]](valX);

Calcular os resíduos

> Residuos:=transform[multiapply[(x,y)-> x-y]]([valY, valYprevisto]);

Gráficos estatísticos - statplots

> ### WARNING: the statplots sub-package has been completely rewritten; see the help pages for details
with(statplots); # certifique: package stats já carregado?

> ValX := [4.535, 4.029, 5.407, 1.605, 5.757, 3.527, 7.890, 8.159, 6.092, 13.442, 2.845, 5.172, 3.277, 8.810, 3.657, 7.226, 3.851, 2.162, 2.668, 4.692];

> ValY := [7.454, 4.476, 2.873, 5.476, 9.975, -1.476, 1.033, 1.140, 4.813, 0.450, -0.788, 9.389, 4.811, -3.107, 4.407, 5.534, 1.691, -0.789, 1.684, 1.605];

> scatterplot( ValX, ValY, color=black );

Boxplot : é um box com

- uma linha central mostrando a mediana,
- uma linha inferior mostrando o primeiro quartil,
- uma linha superior mostrando o terceiro quartil;

> boxplot( ValX, color=red,thickness=2);

>

> with(plots);

Warning, the name changecoords has been redefined

Com display() podemos exibir vários gráficos ao mesmo tempo

> display( { scatterplot( ValX, ValY, color=black),
boxplot( ValY, color=green ), xyexchange(
boxplot( ValX, color=blue , format=notched) ) },
view=[0..17, -4..14], axes=frame );

>

Histogramas

> histogram(ValX, color=green,numbars=6); # soma das áreas=1
# a área de cada barra é a mesma

>

> histogram([random[normald](500)]);
# curva normal aleatória, 500 ptos.

>

> dados2:=
[ -1.96, -.814, 1.86, 1.96, .519, .739, -.0540, .702, .663,
.591, .580, .475, .589, -1.33, .0420, -.460, -.482, 1.58,
.778, .530, -.507, -.233, -.195, .193, -.136];
# dados com distribuição uniforme

Calcular e mostrar FREQÜÊNCIAS

> histogram(dados2, area=count,numbars=6);
# para obter as freqüências, colocar a opção area=count
# a largura de cada barra é a mesma

>

> histogram(dados2, color=yellow);
plot(stats[statevalf,pdf,normald], -3..3, color=red);
plots[display]({%,%%}); # último e penúltimo valores

>

Média, desvio, covariância, com describe()

> with(stats); # certifique: package stats já carregado?
with(describe);

Warning, these names have been redefined: anova, describe, fit, importdata, random, statevalf, statplots, transform

> dados := [1.1, 5.8, 3.4, 4.2, 3.9, 5, 0.9, 6.2];
dados2 := [2.1, 5.8, 4.1, 4.2, 2.9, 5, 1.1, 6.7];
histogram(dados,numbars=8,area=count);boxplot(dados, thickness=3,color=red);

> scatterplot(dados,dados2,color=red,thickness=10);

>

> media:=mean(dados);

> mediana:=median(dados);

> varianc:=variance(dados);

> desviop:=standarddeviation(dados);

> covariance(dados,dados2); # deve haver o mesmo número de valores

Cálculo de densidade probabilística

> with(stats);with(statevalf);

Warning, these names have been redefined: anova, describe, fit, importdata, random, statevalf, statplots, transform

A seguir calcular a

com a média e o desvio padrão calculados antes, (cdf é cumulative density function)

> pr1 := statevalf[cdf,normald[3.8, 1.8] ](2.44);

Agora calcular Pr(2.44<x e x<=10)

>

> pr2 := statevalf[ cdf, normald[3.8, 1.8] ](10) - pr1;

Importação de arquivo de dados

O arquivoX.txt contém os seguintes dados

# após qualquer # é um comentario
# x y
2.1 9.3
3.2 9.0
4.1 * # falta um valor de y, que deve ser indicado com *
5.3 11.0
6.2 21.0
7.7 19.7
8.2 8.1
# fim fim

SE VOCÊ ESTIVER USANDO MS-windows :

Este arquivo deve estar no diretório f:\, onde f:\ seria o nome do diretório.

SE VOCÊ ESTIVER USANDO UNIX:

Este arquivo deve estar no diretório em que este seu worksheet .mws estiver armazenado.

> T:=importdata(`F:\\arquivoX.txt`,2); # note a barra dupla \\ São lidos 2 vetores em T

> scatterplot(T,color=black,thickness=4); # T[1] no eixo x e T[2] no eixo y

> histogram(T[1],color=green,area=count,numbars=6); # só do vetor X

> boxplot(T[1]);

Distribuições de probabilidade: binomial, Poisson, normal, etc.

> ### WARNING: the statplots sub-package has been completely rewritten; see the help pages for details
with(stats[statplots]);

Warning, these names have been redefined: boxplot, histogram, scatterplot, xscale, xshift, xyexchange, xzexchange, yscale, yshift, yzexchange, zscale, zshift

Distribuições:

(1) DISCRETAS:

binomiald[n,p] discreteuniform[a,b]
empirical[list_prob] hypergeometric[N1, N2, n]
negativebinomial[n,p] poisson[mu]

(2) CONTÍNUAS:

beta[nu1, nu2] cauchy[a, b] chisquare[nu]
exponential[alpha, a] fratio[nu1, nu2] gamma[a, b]
laplaced[a, b] logistic[a, b] lognormal[mu, sigma]
normald[mu, sigma] studentst[nu] uniform[a, b]
weibull[a, b]

> y:=[stats[random, binomiald[50,0.6]](20)]; # n=50, p=0.6

> boxplot( y,color=red,thickness=3);

>

> histogram(y,color=red,numbars=10);

>

> y:=[stats[random, normald[50.0,1.0]](20)]; #média 50.0, desvio padrão 1.0

> boxplot( y,color=red,thickness=3);

>

> histogram(y,color=red,numbars=10);

>

> y:=[stats[random, poisson[50]](20)]; #média 50.0

> histogram(y,color=red,numbars=10);

>

> boxplot( y,color=red,thickness=3);

>

The following discrete distributions are available:

binomiald[n,p] discreteuniform[a,b]
empirical[list_prob] hypergeometric[N1, N2, n]
negativebinomial[n,p] poisson[mu]

The following continuous distributions are available:

In the following, the discrete distributions have probability density functions that are evaluated at integral values of x.
The binomiald[n,p] distribution has the probability density function binomial(n,x)*p^x*(1-p)^(n-x). The name binomiald is so chosen in order to distinguish it from the function binomial(). Constraints: x is a non-negative integer no greater than n, n is a positive integer, p is a number between 0 and 1.
The discreteuniform[a,b] distribution has a probability density function that is equal to zero if x<a or if x>b, and equal to 1/(b-a+1) otherwise. Constraints: x is an integer, a<=b
The empirical[list_prob] has its probability density function equal to zero if x<1 or x>nops(list_prob) and equal to list_prob[x] otherwise. Constraints: the probabilities must add to 1 exactly.
The hypergeometric[N1, N2, n], with N1 equal to the size of the success population, N2 equal to the size of the failure population and n equal to the sample size, has the probability density function binomial(N1,x)*binomial(N2,n-x)/binomial(N1+N2,n). Constraints: n<=N1+N2.
The negativebinomial[n, p] distribution has the probability density function equal to

binomial(n+x-1,x)*p^n*(1-p)^x

Constraints: x is a non-negative integer no greater than n, n is a positive integer, p is a number between 0 and 1.

The poisson[mu] distribution has the probability density function exp(-mu)*mu^x/x!
For the continuous distributions, the parameter x takes a real value.
The beta[nu1, nu2] distribution has the probability density function

1/Beta(nu1, nu2) * x^(nu1-1) * (1-x)^(nu2-1).

Constraints: nu1, nu2 are positive integers.

The cauchy[a, b] distribution has the probability density function 1/(Pi*b*(1+((x-a)/b)^2)), b>0.
The chisquare[nu] distribution has the probability density function

x^((nu-2)/2) exp(-x/2)/2^(nu/2)/GAMMA(nu/2), x>0, nu>0.

Constraint: nu is a positive integer.

The exponential[alpha, a] distribution has the probability density function equal to alpha*exp(-alpha*(x-a)) if x>=a and equal to zero if x<a. Constraint: alpha is a positive real number. Default: a=0.
The fratio[nu1, nu2] distribution has the probability density function

GAMMA( (nu1+nu2)/2)/GAMMA(nu1/2)/GAMMA(nu2/2)*(nu1/nu2)^(nu1/2)*
f^((nu1-2)/2) / ( 1+ (nu1/nu2)*f) ^ ((nu1+nu2)/2), x>0, nu1>0, nu2>0

This distribution is also known as the Fisher f distribution and the variance ratio distribution. Constraints: nu1, nu2 are positive integers.

The gamma[a, b] distribution has the probability density function x^(a-1)*exp(-x/b)/GAMMA(a)/b^a, x>0, a>0, b>0. The parameter b, if absent, defaults to the value 1.
The laplaced[a, b] distribution has the probability density function 1/(2*b) * exp(-abs(x-a)/b), b>0. The name laplaced is so chosen to distinguish it from the laplace() function.
The logistic[a, b] distribution has the probability density function

exp( - (x-a)/b ) / b / (1+ exp(-(x-a)/b))^2, b>0

The lognormal[mu, sigma] has the probability density function

1/(x*sqrt(2*Pi)*sigma) * exp( - ((ln(x)-mu)/sigma)^2 /2 ), x>0

The parameter mu has the default value 0 and the parameter sigma has the default value 1. Constraint: sigma cannot be 0. See also the normald distribution.

The normald[mu, sigma] distribution has the probability density function

exp( -(x-mu)^2/2/sigma^2 ) / sqrt(2*Pi*sigma^2)

The parameter mu has the default value 0 and the parameter sigma has the default value 1. Note that sigma is the standard deviation and not the variance. Constraint: sigma must be positive.

The studentst[nu] distribution has the probability density function

GAMMA( (nu+1)/2 )/GAMMA(nu/2)/sqrt(nu*Pi)/(1+t^2/nu)^((nu+1)/2)

Constraint: nu is a positive integer.

The uniform[a, b] distribution has the probability density function equal to 1/(b-a) if a<=x<=b, and to 0 otherwise. The value of b defaults to 1+a. The value of a defaults to 0. Constraint: a<=b.
The weibull[a, b] distribution has the probability density function

a*x^(a-1)*exp( -(x/b)^a )/ b^a, x>0, a>0, b>0

Maple ^TM is a registered trademark of Waterloo Maple Inc.
Math rendered by WebEQ