martes, 29 de mayo de 2012

Octave

Es una herramienta que es parecida a MATLAB, pero es gratis!. Con ella puedes realizar calculos, trabajar con matrices, plotear graficas, entre muchas otras cosas.

Tambien puedes hacer programas con las funciones usuales como: if, while, do entre otras. Puedes salvar tu trabajo realizado o esribir tus instrucciones en un archivo con extension .m, igual que en MATLAB, lo cual los hace compatibles.

Te presento dos ejemplos:

Ejemplo 1.


clear all;
verbose = 1;
mu1 = 4;
mu0 = 1;
N = 2000;
x0 = 0.8;
axis([1 4 0 1])
xlabel "r";
% ylabel "Equilibrium Point, x*";
for i = (1 : N + 1)
  mu = mu0 + (i - 1) * (mu1 - mu0) / N;
  clear x;
  x(1) = x0;
  j = 2;
  notconverged = 1;
  while (notconverged)
    x(j) = mu * x(j - 1) * (1 - x(j - 1));
    j++;
    for k = (1 : j - 2)
      if (abs(x(k) - x(j-1)) < 0.0001)
kstar = k;
notconverged = 0;
break;
      endif
    endfor
  endwhile
  clear muon;
  clear xstar;
  muon = 1;
  xstar = 0;
  for n = (kstar : j - 2)
    muon = [muon ; mu];
    xstar = [xstar ; x(n)];
  endfor
  M = [muon xstar];
  hold on
  gplot M title "" with dots 0
  if (verbose)
    if ((j - 1 - kstar) == 3)
      ans = input("Ready? Enter 1 to continue to the end ...\n");
      if (ans == 1) verbose = 0; endif
    endif
  endif
endfor
%title ("BIFURCATION PLOT FOR THE LOGISTIC MAP (1 < r < 4)")
refresh;
hold off;

* Si escribes estas instrucciones, graficaras una bifurcacion, que se vera como estas:


Ejemplo 2.

% This instructions solve the discrete Larva-Pupa-Adult model of flour beetles.

b=11.677;
Ul=0.5129;
Cea=0.011;
Cel=0.0093;
Cpa=0.0178;
Ua=0.11; %input('input adult mortality Ua:    ')
L0=70;  %input('input initial population L0:     ')
P0=30;
A0=65;
n=40;  %input('end of time interval b:     ')
L=zeros(n+1,1);
P=zeros(n+1,1);
A=zeros(n+1,1);
t=zeros(n+1,1);
L(1)=L0;
P(1)=P0;
A(1)=A0;
for i=1:n
t(i)=i-1;
L(i+1)=b*A(i)*exp(-Cea*A(i)-Cel*L(i));
P(i+1)=L(i)*(1-Ul);
A(i+1)=P(i)*exp(-Cpa*A(i))+A(i)*(1-Ua);
end
t(n+1)=n;
plot(t,L,'-o',t,P,'r:.',t,A,'-*')
legend('Larvae','Pupae','Adults',0)

Y con estas instrucciones veras el crecimeinto de una poblacion adulta de abejas.
Referencias:


lunes, 23 de abril de 2012

R 64


R es un lenguaje y ambiente para estadistica computacional y gráficas. Es un proyecto GNUsimilar al ambiente y lenguaje S, el cual fue desarrollado en los Laboratorios Bell (anteriormente AT&T, ahora “Lucent Technologies”) por John Chambers y sus colegas. R puede ser considerado como una implementación diferente de S. Hay algunas diferencias importantes, pero la mayoría del código escrito para S corre en R.

R provée una amplia variedad de estadísticas (modalción lineal y nonlineal, pruebas clásicas de estadística, análisis de series de tiempo, clasificación, agrupamiento, ...) y técnicas gráficas, y es altamente extensible. El lenguaje S es amenudo un vehículo a escoger para investigación en metodologías estadísticas, y R es “Open Source”.

Una de la fortalezas de R es la facilidad de realizar gráficas bien diseñadas con calidad de publicación, incluyendo símbolos matemáticos y fórmulas donde se necesiten.


Ejercicios:



x <- rnorm(10)
y <- rnorm(10)
plot(x,y)
plot(x,y, xlab="Diez numeros al azar", ylab="Otros diez numeros", xlim=c(-2, 2), ylim=c(-2,2), pch=22, col="red", bg="yellow", bty="l", tcl=0.4, main="Como personalizar un grafico en R", las=1, cex=1.5)

opar <- par() ##abre una ventana en R64
par(bg="lightyellow", col.axis="blue", mar=c(4,4,2.5,0.25))
plot(x,y, xlab="Diez numeros al azar", ylab="Otros diez numeros", xlim=c(-2, 2), ylim=c(-2,2), pch=22, col="red", bg="yellow", bty="l", tcl=-0.25, las=1, cex=1.5) ## pone las figuras, pone ejes, escalas y dibuja fondo
title("Como personalizar un grafico en R (bis)", font.main=3, adj=1) ## pone titulo.


##se puede sobreescribir las 2 ultimas lineas y obtenemos fig5b
par(bg="lightgrey", col.axis="blue", mar=c(4,4,2.5,0.25))
plot(x,y, xlab="Diez numeros al azar", ylab="Otros diez numeros", xlim=c(-2, 2), ylim=c(-2,2), pch=22, col="cyan", bg="pink", bty="l", tcl=-0.25, las=1, cex=1.5)

opar <- par()
par(bg="lightgrey", mar=c(2.5, 1.5, 2.5, 0.25))
plot(x,y, type="n", xlab="", ylab="", xlim=c(-2, 2), ylim=c(-2,2), xaxt="n", yaxt="n")
rect(-3, -3, 3, 3, col="cornsilk")
points(x, y, pch=10, col="red", cex=2)
axis(side=1, c(-2, 0, 2), tcl=-0.2, labels=FALSE)
axis(side=2, -1:1, tcl=-0.2, labels=FALSE)
title("Como personalizar un grafico en R (ter)", font.main=4, adj=1, cex.main=1)
mtext("Diez números al azar", side=1, line=1, at=1, cex=0.9, font=3)
mtext("Otros diez números", line=0.5, at=-1.8, cex=0.9, font=3)
mtext(c(-2, 0, 2), side=1, las =1, at=c(-2,0,2), line=0.3, col = "blue", cex=0.9)
mtext(-1:1, side=2, las=1, at=-1:1, line=0.2, col = "blue",  cex=0.9)


LATTICE y GRID

LATTICE Y GRID…
library(lattice)
n <- seq (5, 45, 5)
> n
[1]  5 10 15 20 25 30 35 40 45
y <- factor(rep (n,n), labels=paste("n =", n))
> y
  [1] n = 5  n = 5  n = 5  n = 5  n = 5  n = 10 n = 10 n = 10 n = 10 n = 10 n = 10
 [12] n = 10 n = 10 n = 10 n = 10 n = 15 n = 15 n = 15 n = 15 n = 15 n = 15 n = 15
 [23] n = 15 n = 15 n = 15 n = 15 n = 15 n = 15 n = 15 n = 15 n = 20 n = 20 n = 20
 [34] n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 20
 [45] n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 25 n = 25 n = 25 n = 25 n = 25
 [56] n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25
 [67] n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 30 n = 30
 [78] n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30
 [89] n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30
[100] n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 35 n = 35 n = 35 n = 35 n = 35
[111] n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35
[122] n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35
[133] n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 40 n = 40 n = 40
[144] n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40
[155] n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40
[166] n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40
[177] n = 40 n = 40 n = 40 n = 40 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45
[188] n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45
[199] n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45
[210] n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45
[221] n = 45 n = 45 n = 45 n = 45 n = 45
Levels: n = 5 n = 10 n = 15 n = 20 n = 25 n = 30 n = 35 n = 40 n = 45
densityplot(~x|y)


data(quakes)
mini <- min(quakes$depth)
maxi <- max(quakes$depth)
int <- ceiling((maxi-mini)/9)
inf <- seq(mini, maxi, int)
quakes$depth.cat <- factor(floor(((quakes$depth - mini) / int)), labels=paste(inf, inf + int, sep="-"))
xyplot(lat ~ long | depth.cat, data = quakes)


Lee los datos que son pruebas de 6 diferentes tipos de insecticidas, utilizando el numero de insectos muertos como la variable de respuesta. Cada insecticida se probo 12 veces, para un total de 72 observaciones Con esos datos realizaremos un analisis de varianza simple de respuesta como funcion del insecticida usado. Primero se deben cargar los datos en la memoria (con data) y el analisis se realiza con la funcion aov().

data(InsectSprays)
aov.spray <- aov(sqrt(count) ~ spray, data = InsectSprays)

El argumento principal (y obligado) de aov() es una formula que especifica la respuesta en el lado izquierdo del símbolo ~ hy la variable explicativa en el lado derecho. La opción data = InsectSprays indica que las variables deben ser buscadas en el marco de datos InsectSprays. Est sintaxis es equivalente a:

aov.spray <- aov(sqrt(InsectSprays$count) ~ InsectSprays$spray)
o también, si conocemos el numero de las columnas de las variables:
aov.spray <- aov(sqrt(InsectSprays[, 1]) ~ InsectSprays[,2])

Los resultados del análisis no se muestran ya que son asignada a un objeto llamado aov.spray. Usaremos otras funciones para extraer los resultados, como print() para ver un resumen breve del análisis (sobre todo los parámetros) y sumary() para ver mas detalles (incluyendo pruebas estadísticas):

aov.spray
Call:
   aov(formula = sqrt(count) ~ spray, data = InsectSprays)

Terms:
                   spray Residuals
Sum of Squares  88.43787  26.05798
Deg. of Freedom        5        66

Residual standard error: 0.6283453 
Estimated effects may be unbalanced

summary(aov.spray)
            Df Sum Sq Mean Sq F value Pr(>F)    
spray        5  88.44  17.688    44.8 <2e-16 ***
Residuals   66  26.06   0.395                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Recordemos que escribir el nombre del objeto como un comando, es similar al comando print(aov.spray). Para ver una representación gráfica de los resultados podemos usar plot() o termplot(). Antes de escribir plot(aov.spray) dividiremos la ventana gráfica en cuatro partes de tal manera que las cuatro gráficas diagnosticas se dibujan en la misma ventana. Los comandos son:
> opar <- par()
> par(mfcol = c(2,2))
> plot(aov.spray)
y se obtiene la grafica 13.
y con la instrucción
termplot(aov.spray, se=TRUE, partial.resid=TRUE, rug=TRUE)se obtiene el grafico 14 

Funciones como aov o lm devuelven una lista donde los diferentes elementos corresponden a los resultados del análisis. Podemos ver la estructura del objeto devuelto por aov() con la instrucción:
str(aov.spray, max.level = -1)
List of 13
 - attr(*, "class")= chr [1:2] "aov" "lm"

Otra manera de ver esta estructura es visualizando los nombres del objeto, con la instrucción:
names(aov.spray)
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "contrasts"     "xlevels"       "call"          "terms"        
[13] "model"   

Los elementos se pueden exraer como en cualquier otra lista, con la instrucción:
aov.spray$coefficients
(Intercept)      sprayB      sprayC      sprayD      sprayE      sprayF 
  3.7606784   0.1159530  -2.5158217  -1.5963245  -1.9512174   0.2579388 

summry() también crea una lista la cual, en el caso de aov(), es simplemente una tabla de pruebas:
str(summary(aov.spray))
List of 1
 $ :Classes ‘anova’ and 'data.frame': 2 obs. of  5 variables:
  ..$ Df     : num [1:2] 5 66
  ..$ Sum Sq : num [1:2] 88.4 26.1
  ..$ Mean Sq: num [1:2] 17.688 0.395
  ..$ F value: num [1:2] 44.8 NA
  ..$ Pr(>F) : num [1:2] 6.33e-20 NA
 - attr(*, "class")= chr [1:2] "summary.aov" "listof"
> names(summary(aov.spray))
NULL


Las funciones genéricas también se denominan métodos. Esquematicamente, se construyen como methodfoo, donde foo es la función de análisis. En el caso de summary, podemos ver las funciones que se usan en este método.
apropos("^summary")
 [1] "summary"                 "Summary"                 "summary.aov"            
 [4] "summary.aovlist"         "summary.connection"      "summary.data.frame"     
 [7] "Summary.data.frame"      "summary.Date"            "Summary.Date"           
[10] "summary.default"         "Summary.difftime"        "summary.factor"         
[13] "Summary.factor"          "summary.glm"             "summary.infl"           
[16] "summary.lm"              "summary.manova"          "summary.matrix"         
[19] "summary.mlm"             "Summary.numeric_version" "Summary.ordered"        
[22] "summary.POSIXct"         "Summary.POSIXct"         "summary.POSIXlt"        
[25] "Summary.POSIXlt"         "summary.srcfile"         "summary.srcref"         
[28] "summary.stepfun"         "summary.table"           "summaryRprof"    

Podemos ver la diferencia para este método comparando como actúa en una regresión lineal contra un análisis de varianza:
x <- y <- rnorm(5);
mod <- lm(y ~ x)
names(mod)
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model" 
names(summary(mod))
 [1] "call"          "terms"         "residuals"     "coefficients" 
 [5] "aliased"       "sigma"         "df"            "r.squared"    
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled" 

Los objetos devueltos por aov(), lm(), summary(), … son listas, pero no se visualizan como las listas "comunes y corrientes" que hemos visto anteriormente. De hecho, son métodos sprint de estos objetos (recuerde que escribir el nombre de un objeto equivale a usar print()):
apropos("^print")
 [1] "print"                       "print.anova"                
 [3] "print.AsIs"                  "print.by"                   
 [5] "print.coefmat"               "print.condition"            
 [7] "print.connection"            "print.data.frame"           
 [9] "print.Date"                  "print.default"              
[11] "print.density"               "print.difftime"             
[13] "print.DLLInfo"               "print.DLLInfoList"          
[15] "print.DLLRegisteredRoutines" "print.factor"               
[17] "print.family"                "print.formula"              
[19] "print.ftable"                "print.function"             
[21] "print.glm"                   "print.hexmode"              
[23] "print.hsearch"               "print.infl"                 
[25] "print.integrate"             "print.libraryIQR"           
[27] "print.listof"                "print.lm"                   
[29] "print.logLik"                "print.NativeRoutineList"    
[31] "print.noquote"               "print.numeric_version"      
[33] "print.octmode"               "print.packageInfo"          
[35] "print.POSIXct"               "print.POSIXlt"              
[37] "print.proc_time"             "print.restart"              
[39] "print.rle"                   "print.simple.list"          
[41] "print.srcfile"               "print.srcref"               
[43] "print.summary.table"         "print.summaryDefault"       
[45] "print.table"                 "print.terms"                
[47] "print.ts"                    "print.warnings"             
[49] "printCoefmat"               
Todos estos métodos de print permiten cierta visualización dependiendo del análisis.

Funcion que realiza el análisis numérico del modelo de Ricker. Este modelo se usa ampliamente en ecología, particularmente en el estudio demográfico de peces. El objetivo es simular el modelo con respecto a la tasa de crecimiento r y el numero inicial de individuos en la población N_0 (la capacidad de carga K es fijada en 1), los resultados se mostraran como una gráfica del numero de individuos en funcion del tiempo. Agregaremos una opción ver los resultados en los últimos pasos de la simulación (por default los resultados son graficados).
> ricker <- function(nzero, r, K=1, tiempo=100, desde=0, hasta=tiempo)
+ { N<- numeric(tiempo+1)
+ N[1] <- nzero
+ for(i in 1:tiempo) N[i+1] <- N[i]*exp(r*(1-N[i]/K))
+ Tiempo <- 0:tiempo
+ plot(Tiempo, N, type="l", xlim = c(desde, hasta))
+ }
> layout(matrix(1:3, 3, 1))
> ricker(0.1, 1); title("r=1")
> ricker(0.1, 2); title("r=2")
> ricker(0.1, 3); title("r=2")




Bibliografia:



martes, 27 de marzo de 2012

Phyton V2.7

Se trata de comparar diferentes formas de hacer lo mismo.

1. Hacer un programa en lenguaje C, que lea un archivo de entrada y busque un string y si lo encuentra, lo guardara en una variable con valor 1.

/ * Programa que encuenta una secuencia de letras de en un archivo que abre *
">*=================================================*

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int main(int argc, char *argv[])
{
int cont=0;

char renglon [50];
FILE *ifp;
if (argc != 2) {"\n%s%s%s\n\n", "Debes teclear: ", argv[0], " ArchivoEntrada");"r");while ((fgets (renglon , 50 , ifp)) != NULL){
if(strstr(renglon,"NORMAL COMPLETION") != NULL)
cont++;
}
printf (
"\n\n%d\n\n", cont);

return cont;
} 
Se compila el programa y no tiene errores, entonces se hace un archivo ejecutable para que lo repita n veces (10 en mi caso)
archivo EnPal.sh





#!/bin/bash
echo "Buscara una cadena de caracteres en 10 archivos."
for j in 20; do 
for i in 2; do
for k in 1 2 3 4 5; do
# echo "wee_"$j"_"$i"_"$k
./ex "wee_"$j"_"$i"_"$k".lst" >> "ResEncPal.txt"
done
done
#awk '{if ($1 == 1)} {print "El archivo wee_"$j"_"$i"_"$k".lst fue exitosos"}' "ResEncPal.txt"
done 
#awk '{if ($1 == 1)} {print "El archivo wee_"$j"_"$i"_"$k".lst fue exitosos"}' "ResEncPal.txt"
for j in 40; do 
for i in 4; do
for k in 1 2 3 4 5; do
# echo "wee_"$j"_"$i"_"$k
./ex "wee_"$j"_"$i"_"$k".lst" >> "ResEncPal.txt"
done
done
done 
echo "F I N" 

El archivo ResEsnPal.txt, solo contiene 0 y 1's, 1: si el archivo que lee contiene la palabra NORMAL y 0 en otro caso.

2. Con Python.

Primero concatenamos los 10 archivos con la instruccion cat (cat arch1 arch2 ... arch10 > prov) y despues ejecutamos las siguientes lineas: 

python
mport string
text = open('prov').read()
print string.count(text, 'NORMAL')
quit()



3.  Desde Unix 
 
 
grep -c 'NORMAL' prov



N -> es el numero de veces que el status fue terminacion completa


Referencias:

Referencia 1.

Gnuplot es una herramienta muy versatil para hacer graficas en 2 y 3 dimensiones.

A continuacion algunos ejemplos:


Genera 2 curvas y una mas con linea rellenas:set title
set key outside
plot [-10:10] [-5:3] 1.5+sin(x)/x with filledcurves, sin(x)/x with filledcurves, 1+sin(x)/x with lines


Genera 2 parabolas encontradas:
set key on
set title "Intersection of 2 parabolas"
plot x*x with filledcurves, 50-x*x with filledcurves, x*x with line lt 1
Genera el seno y coseno rellenos, con un grid:set title "Filled sinus ans cosinus curves"
set grid front
set title "Filled sinus ans cosinus curves"
plot 2+sin(x)**2 with filledcurve x1, cos(x)**2 with filledcurve x1
Genera una red de 1000 numeros aleatorios:unset key
set xrange [0:1]
set yrange [0:1]
set zrange [0:1]
set title "Lattice test for random numbers"
set xlabel "rand(n) ->"
set ylabel "rand(n + 1) ->"
set zlabel "rand(n + 2) ->"
set format x "3.2f"
set format y "3.2f"
set format z "3.2f"
set tics
set sample 1000
set style function impulses
set parametric
splot rand(0), rand(0), rand(0)
Grafica en 3 D una malla con diferentes colores... set title "Textcolor options" tc lt 1
set samples 20; set isotemplates 20
set autoscale
set key box
set pm3d at s
set color box horizon user origin .1, .12, size .8, .015
set view 58, 64, 0.83
set xrange [-10:10]
set yrange [-10:10]
set zrange [-10:10]
#Test lables
set label 1 "textcolor palette z" at 12, -10, -10 nopoint tc pal z
set label 3 "tc pal z" at 12, -6, -6 nopoint tc pal z
set label 4 "tc pal z" at 12, -3, -3 nopoint tc pal z
set label 10 "textcoor lt 1" at -10, -8, 24 front nopoint tc lt 1
set label 10 "textcoor lt 1" at -10, -8, 24 front nopoint tc lt 1
set label 5 "tc pal z" at 12, -0, 0  nopoint tc pal z
set label 6 "tc pal z" at 12, 3 , 3  nopoint tc pal z
set label 7 "tc pal z" at 12,  6,  6 nopoint tc pal z
set label 8 "tc pal z" at 12, 9 , 9  nopoint tc pal z
#
set label 'xlabel should be lt 4' tc lt 4
set cblabel 'color cblabel' textcolor lt 3
#
set label 10 "textcolor lt 1" at -10, -8, 24 front nopoint tc lt 1
set label 11 "textcolor lt 2" at -10, -8, 21 front nopoint tc lt 2
set label 12 "textcolor lt 3" at -10, -8, 18 front nopoint tc lt 3
set label 13 "textcolor default" at -10, -8, 15 front nopoint tc def
set label 14 "textcolor cb 5" at -10, -8, 12 front nopoint tc pal cb 5
set label 15 "tc cb 0" at -10, -8, 9 front nopoint tc pal cb 0
set label 16 "tc cb -5" at -10, -8, 6 front nopoint tc pal cb -5
set label 17 "textcolor frac .75" at -10, -8, 3 nopoint tc pal frac .75
set label 18 "tc frac .25" at -10, -8, 0 nopoint tc pal frac .25
#
set ylabel 'ylabel should still be black'
splot y

Para mandar el dibujo a un archivo:

a)set terminal png        # gnuplot recomends setting terminal
                            before output
set output "output.png"   # the output filename; to be set after
                            setting terminal
replot                    # es para dibujar sobre una grafica

Manda el archivo a una hoja de presentacion... salir de gnuplot antes de buscarlo y cambiar a directorio diana
set term post enhanced
set out 'gplt.eps'
set xlabel '{/symbol q_1})'
set ylabel 'sin^2({/symbol q_1})
plot sin(x)**2

Imprimime la funcion seno y coseno en la misma gráfica con dos tipos diferentes...
plot sin(x) with linespoints pointtype 5, cos(x) w boxes lt 4Imprime 4 graficas en un solo archivo, dividiendo en 4 partes iguales:set xrange [-pi:pi]
set size 1,1
set origin 0,0
set multiplot
set size 0.5,0.5
set origin 0,0.5
plot sin(x)
set size 0.5,0.5
set origin 0,0
plot 1/sin(x)
set size 0.5,0.5
set origin 0.5,0.5
plot cos(x)
set origin 0.5,0.5
set size 0.5,0.5
set origin 0.5,0
plot 1/cos(x)
unset multiplot

reset


Leer desde archivo.
Se tiene el siguiente archivo llamado, ibm.dat
10-jun-04   90.23 90.75 89.89 90.46
9-jun-04    89.90 90.55 89.81 90.09
8-jun-04    88.64 90.50 88.40 90.04
7-jun-04    88.75 88.99 88.01 88.64
4-jun-04    87.95 88.49 87.50 87.56

3-jun-04    87.85 88.10 87.35 87.35
2-jun-04    88.64 88.64 87.89 87.98
1-jun-04    88.00 88.48 87.30 88.12
#date       open  high  low   close
para graficar el precio de apertura (col. 2) usando el linespoints style para unir con un a linea que interpolara...set xdata time
set timefmt "%d-%b-%y"
set format x "%b %d"
plot ["31-May-04":"11-Jun-04"] 'ibm.dat' using 1:2 with linespoints
para imprimir todos los datos:set xdata time
set timefmt "%d-%b-%y"
set format x "%b %d"
set bars 5
plot ["31-May-04":"11-Jun-04"] 'ibm.dat' using 1:2:3:4:5 with financebars

Para graficar mis datos:
set xrange [0:20]                                                               set yrange[0:30]                                                               set tics scale .5
set xtics ("20.2.1" 1, "20.2.2" 2, "20.2.3" 3, "20.2.4" 4, "20.2.5" 5, "20.3.1" 6, "20.3.2" 7, "20.3.3" 8, "20.3.4" 9, "20.3.5" 10, "40.3.1" 11, "40.3.2" 12, "40.3.3" 13, "40.3.4" 14, "40.3.5" 15, "40.4.1" 16,"40.4.2" 17,"40.4.3" 18,"40.4.4" 19," 40.4.5" 20)
set terminal png (puede ser eps)
set output "grafica_diana.png"
plot 'datos_diana.dat' using 1:3:4:5 title 'Resultados de 20 instancias' with yerrorbars
replot 'datos_diana.dat' using 1:3 title 'valor de la funcion objetivo' w  points 1    
equivalen las ultimas dos lineas a:         
plot 'datos_diana.dat' using 1:3:4:5 title 'Resultados de 20 instancias' with yerrorbars, 'datos_diana.dat' using 1:3 title 'valor de la funcion objetivo' w  points 1
quit

Referencias:


martes, 13 de marzo de 2012

AWK y SED

Se muestra el uso de dos herramientas nuevas:
1. AWK
2. SED

En ambos casos se crearon programas del tipo BASH para realizar repeticiones de los programas con ambas herramientas.

El objetivo de ambos programas es obtener información de un archivo *.lst generado al resolver un programa de optimización en GAMS, cada uno con su sintaxis respectiva.

El primer archivo se llama ultimo.sh y usa AWK par leeer por columnas la informacion del archivo de entrada (*.lst) y cuando encuentra las coincidencias re-direcciona la infornación hacia diferentes archivos de salida (res_weee_NumNodos_Param_Repetición.txt) El segundo archivo se llama intento_sed.sh y usa

El archivo ultimo.sh que pide el numero de nodos para saber que instancias son las que se realizaran (20 o 40 nodos) cada uno con 10 variantes.

El archivo intento_sed.sh pide el número de nodos para crear el nombre del archivo de entrada (*.lst) y cuando encuentra las coincidencias re-direcciona la información hacia diferentes archivos de salida (obtenido_NumNodos_Param_Repetición)

Se adjuntan ambos programas:

El archivo ultimo.sh


#!/bin/bash
echo "numero de unidades basicas: "
read
v=$REPLY
echo $v
for i in  1 2 3 4 5; do
    if [ $v -eq "20" ]
        then
        for j in 2 3; do
awk 'BEGIN {print "R E S U L T A D O S \n"} $1 == "****" && $3 == "STATUS" {print $2,$3,":",$5,$6,$7} $2 == "OBJECTIVE" {print $2,$3,$4} $1 == "MIP" && $2 == "status" {print $0} $1 == "Best" {print $0, "\nLos GAP absolutos y relativos son:"} $2 ~ /gap:$/ {print $0} END {print "--------------fin----------------"}' $"weee_"$v"_"$j"_"$i"_cplex_opt.lst" > "res_weee_"$v"_"$j"_"$i"_cplex_opt.txt"
        done
    elif [ $v -eq "40" ];
        then
        for j in 3 4; do
awk 'BEGIN {print "R E S U L T A D O S \n"} $1 == "****" && $3 == "STATUS" {print $2,$3,":",$5,$6,$7} $2 == "OBJECTIVE" {print $2,$3,$4} $1 == "MIP" && $2 == "status" {print $0} $1 == "Best" {print $0, "\nLos GAP absolutos y relativos son:"} $2 ~ /gap:$/ {print $0} END {print "--------------fin----------------"}' $"weee_"$v"_"$j"_"$i"_cplex_opt.lst" > "res_weee_"$v"_"$j"_"$i"_cplex_opt.txt"
        done

    else
        echo "Solo hay instancias con 20 y 40 nodos"
        break
    fi
    echo "Finalizado con $v unidades basicas"
done

El formato del archivo de salida es el siguiente:
R E S U L T A D O S
SOLVER STATUS : NORMAL COMPLETION
MODEL STATUS : OPTIMAL
OBJECTIVE VALUE 2.0000
Best possible:           2.000000
Los GAP absolutos y relativos son:
Absolute gap:            0.000000
Relative gap:            0.000000
--------------fin----------------


El archivo intento_sed.sh

#!/bin/bash
echo "numero de unidades basicas: "
read
v=$REPLY
for i in  1 2 3 4 5; do
    if [ $v -eq "20" ];
        then
        for j in 2 3; do
            sed -n -e '/STATUS/p' -e '/gap/p' "weee_"$v"_"$j"_"$i"_cplex_opt.lst" > "obtenido_"$v"_"$j"_"$i
        done
    elif [ $v -eq "40" ];
        then
        for j in 3 4; do
            sed -n -e '/STATUS/p' -e '/gap/p' "weee_"$v"_"$j"_"$i"_cplex_opt.lst" > "obtenido_"$v"_"$j"_"$i
        done
    else
        echo "Solo hay instancias con 20 y 40 nodos"
        break
    fi
    echo "Finalizado con $v unidades basicas"
done

El formato del archivo de salida es el siguiente:
**** SOLVER STATUS     1 NORMAL COMPLETION
**** MODEL STATUS      1 OPTIMAL
>>  epgap  0.0001
Absolute gap:            0.000000
Relative gap:            0.000000