viernes, 6 de enero de 2012

interpolacion wavelet


Un estudio del Escalamiento de Imágenes mediante Interpolación basada en la Transformada Discreta Wavelet de Daubechies

La interpolación es una poderosa herramienta matemática cuya versatilidad ha permitido que sea aplicada sobre imágenes. Entre las principales áreas a las que se aplica la interpolación de imágenes son: computación gráfica, renderizado, edición de imágenes,  reconstrucción de imágenes médicas, visualización de imágenes online, entre otras.  Interpolación trabaja con datos conocidos para así aproximar datos desconocidos por ejemplo:  al interpolar los datos de color azul obtenemos una función aproximada sin(x) e interpolamos para obtener nuevos puntos (rojo )
Interpolación de la función Seno
Interpolación de imágenes consiste en redimensionar la imagen ya sea en zooming o en scaling . Por ejemplo en la siguiente Figura
Interpolación 2-D de una imagén 

Las técnicas de interpolación de imágenes se agrupan en dos categorías: técnicas adaptativas y no adaptativas. Los principios de los algoritmos de interpolación adaptativa se basan en las características intrínsecas de la imagen a diferencia de los algoritmos no adaptativos que no consideran esto y repiten un procedimiento por cada píxel o grupos de píxeles de una imagen. 

Siendo la Transformada Discreta de Wavelets (DWT) una técnica utilizada en diferentes aplicaciones de procesamiento digital de imágenes debido a que brinda mejores herramientas para realizar análisis multi-resolución de señales y basados en los trabajos de [Acharya, 2006] y [Acharya, 2007], en este documento presentamos los resultados obtenidos al realizar el escalamiento de imágenes en tres diferentes espacios de color (RGB, CIE Lab y HSV) utilizando interpolación basada en la DWT de Daubechies de orden 4.

Una transformada wavelet


y su transformada inversa

La transformada discreta de wavelets (DWT) se basa en el uso de valores discretos para los parámetros de escala s y traslación τ. La DWT usa bancos de filtros para la construcción de un plano multi-resolución en el dominio de tiempo y frecuencia. Para el análisis y reconstrucción de señales, la DWT utiliza bancos de filtros multi-resolución y un filtro de wavelet especial.

Bancos de filtros consisten en filtros que separan una señal en bandas de frecuencias, una señal x(k) es filtrada por L(z) (filtro paso bajo) y H(z) (filtro paso alto).  

Donwsampling: Construccion de la señal en un factor de 1/2 ( ↓2)
Upsampling: Construcción de la señal en un factor de 1/2 (↑ 2)

banco de filtro de dos canales

Banco de filtros en tres niveles
 Interpolación Wavelet

 Es una aplicación de la transformada wavelet con un filtro en especial, para nuestro caso usamos el filtro daubechies 4 .
Descomposición en 3 niveles con 2D DWT. (a) Descomposición en el primer nivel. (b) Descomposición en el segundo nivel. (c) Descomposición en el tercer nivel
Como los filtros waveletes son separables DWT 2-D se pueden implementar en DWT 1-D aplicando a las filas y despues a las columnas o a la inversa. En la Figura anterior es la descomposicion en 3 niveles de una imagen usando filtros de paso bajo y alto.



La interpolación wavelet propuesta por [Acharya, 2007] es definida en la siguiente imagen



Proceso de upsampling para una imagen. (a) Forward DWT de la imagen de entrada. (b) construcción de la imagen DWT virtual. (c) Imagen escalada a una ampliación después de aplicar la DWT inversa.


Consta de la siguientes etapas:
  1. Descomposición de la imagen I en un nivel asi obtenemos las bandas LL, HL, LH ,HH
  2. Creación de la imagen I' con el doble de tamaño de la imagen original (2mx2n)
  3. Logicamente se divide I' en 4 cuadrantes en el que:
    • El primer cuadrante es la imagen I multiplicada por un escalar s , segun [Acharya,2007] s es ganacia al cuadrado del filtro, pero para nuestro caso decidimos que sea un factor escalar variante y no dependiente.
    • En el segundo y tercer cuadrante se realiza el completado con cero de los pixeles intermedios
    • El último  cuadrante su valor es cero.
  4. Aplicamos la transformada wavelet inversa y obtenemos nuestra imagen inteporlada :D .
Para nuestro experimento usamos una base de imagenes con 2500 muestras obtenidas de http://staff.science.uva.nl/~aloi/ en diferentes modelos de color RGB, CIELab y HSV .  Obteniendo los siguientes resultados:



El siguiente gráfico es un resumén de nuestro experimento


Más detalles en [1]
[1] Calderón Niquín Marks, Valverde-Rebaza Jorge (2011), Un estudio del Escalamiento de Imágenes mediante Interpolación basada en la Transformada Discreta Wavelet de Daubechies, Proceeding pag 85-88 de III Congreso de Computación y Telecomunicaciones - Universidad Inca Garcilazo de la Vega 
Puedes descargarlo en el sgte. link

lunes, 26 de diciembre de 2011

Writing for Computer Science (Springer)

Writing for Computer Science - J. Zobel Springer 1997) WW

Writing for computer science students

Writing for Computer Science students

viernes, 16 de diciembre de 2011

Descargar Geo Dataset desde R

Cada Geo Dataset tiene un codigo o identificador unico por ejemplo en esta web puedes navegar el tipo de datos que deseas o de enfermeda,estudio deseas descargar http://www.ncbi.nlm.nih.gov/gds en nuestro caso buscamos Breast cancer relapse free survival y seleccionamos el item con codigo : GSE2034 http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE2034 que es el caso de estudio de metástasis de cancer mamario.

Bueno en R tienes que hacer los sgtes pasos:

1. Instalamos la libreria "GEOquery"
biocLite("GEOquery")

pero nos aseguramos que este ya instalado el source de bioconductor

2.referenciamos la libreria
library("GEOquery")

3. Descargamos nuestra data
getGEOSuppFiles('GSE2034')

4. si deseas descomprimirlo desde R :
untar("GSE2034_RAW.tar", exdir="data")

Instalar bioconductor en R

Para instalar bioconductor debes seguir los sgtes pasos.

1. Ejecuta a R como administrador si es en linux (sudo R).
2. Añadimos el script de bioconductor


> source("http://www.bioconductor.org/biocLite.R")


3. Llamamos a la sgte funcion del script instalado

> biocLite()

4. Esperamos a que compile :D

5. Verficamos las nuevas librerias instaladas con

> library()

Instalar paquetes en R

Si estas en linux o windows inicia como administrador a R . Por ejemplo en un linux:

> sudo R

una vez abierto R con permisos de superusuario

> install.packages("nombrepaquete",dependencies=true)

dependencies: es true para q instale automatico otros paquetes del cual depende

por ejemplo queremos instalar lo sgte :

> install.packages("graphics",dependencies=TRUE)

y si queremos instalar una lista de paquetes

install.packages(c("MASS","nnet","survival","splines"),dependencies=TRUE)

domingo, 6 de noviembre de 2011

Mi último ACM-ICPC

Mi último ACM-ICPC

EL INICIO

Hace ya 2 años atras o más pensé porque si somos una carrera de Ciencia de la Computación porque no participamos en el ACM-ICPC , las respuestas eran obvias:
  • Poco interes de personas en ello .
  • Poco apoyo para viajar a Lima a participar.
  • Entre otras variables más.
Asi q el 2010 me arriesgue con todo a traer una sede a Trujillo para que mi querida Universidad Nacional de Trujillo en especial mi escuela participe :D . Si que fue un trabajo duro en la UNT tanto en la organizacion del concurso del concurso como buscar personas que participen en el evento. Asi que me meti de asistente de cátedra de un curso de Tecnicas de programación ("diseño de algoritmos" en mi U) y deje de llevar un curso =/ ; en los labs les dejaba ejercicios tipo para el concurso .


De esa experiencia rescate un par de personas y con otros amigos del ultimo año formamos el team ACM-ICPC UNT :D .

Como todo inicio siempre es dificil: sin coach, ni apoyo de autoridades ni profes y nadie que un peso x nosotros =( .

Bueno en ese año 2010 todo nos salio con las justas pq también teniamos q organizar un congreso y las mismas personas q estábamos en los teams eramos organizadores del congreso =/ .

Se llevo a cabo todo :D , yo con un stress pro el cual me fui de viaje x 2 semanas porq ya no daba para más .

Este año

Después de nuestra primera experiencia el 2010, los sobrevivientes de aquel entonces eramos 3 a 4 personas , nadie más =/ . Llego un tiempo q solo mi team quedo , por esos momentos pensé esto se cierra =( no da para mas ='( . Pero gracias a Dios x hay se me aparecio un nuevo team , jovén y con proyección a futuro :D .

Así que en unos 3 a 4 meses entrenamos :D , sin coach x un tiempo hasta recibir la ayuda de mi amigo Nils como coach .

Y bueno no entranamos al nivel q deseamos pero ya algo más q el año pasado :D . No tuvimos financiamento de nuestra universidad para remate =/ , en fin pequeños detalles q nos hacen más fuertes :D .

El concurso

Este año lo organizo al Universidad Privada del Norte , algo a destiempo pero se llevo a cabo , con algún u otro inconveniente pero salio :D . La UNT nos presentamos con dos equipos : UNTrainning (mi team) y Time Limit Exced , bueno solo para dos teams nos alcanzo nuestro dinero .

Este año me pareció que el set de problemas estuvo mas complejo q el del 2010 , intentamos 5 de los cuales solo nos acepto dos xD (ya algo es algo ). Y ganamos localmente con aquellos resueltos . Esperaba más pero x A o B motivos no se pudo. Espero que los teams sucesores continuen con la chamba y algún día estemos en la cima :D .

Ocurrencias del ACM-ICPC

Por este concurso me:
  • Estresee
  • Me endeude el 2010 con 900 lukas :| hasta q me dio el dinero la unt xD .
  • Puse n veces de mi bolsillo para q esto corra.
  • Di el tiempo que no tenia pero = lo daba.
  • Deje de llevar un curso.
  • Hasta con la mafia de la unt(autoridades) xD me involucre para financiamiento en el 2010.

Todos pensarán que webon es este q desperdicia tanto de si . Me da = lo que digan xD , me siento satisfecho de ello y lo volvería a realizar y más :D .

Este año ya no representaré a mi querida universidad ( "entre lagrimas :'( " ) , pero di el segundo paso de los 6 q planee para mi U :D , espero que los que se quedan continuen con este trabajo en vez q me agradezcan :D . También espero el sgte año ganarme el puesto de coach , seguiré entrenando :D weno ya no tendré otras ocupaciones como este año .

Mil gracias a Dios por guiarme en esto :D y a todos por su apoyo :D . Consejo : Si piensas esperar que te den todo o te caiga del cielo , estás perdido y vivirás mal a lo largo de tu vida = / . Hasta que tus sueños se conviertan en metas :D .

Por aca unas fotos :
























viernes, 30 de septiembre de 2011

Introductory Statistics with R


Introductory Statistics with R

This book provides an elementary-level introduction to R, targeting both non-statistician scientists in various fields and students of statistics. The main mode of presentation is via code examples with liberal commenting of the code and the output, from the computational as well as the statistical viewpoint. A supplementary R package can be downloaded and contains the data sets.

The statistical methodology includes statistical standard distributions, one- and two-sample tests with continuous data, regression analysis, one- and two-way analysis of variance, regression analysis, analysis of tabular data, and sample size calculations. In addition, the last six chapters contain introductions to multiple linear regression analysis, linear models in general, logistic regression, survival analysis, Poisson regression, and nonlinear regression.

In the second edition, the text and code have been updated to R version 2.6.2. The last two methodological chapters are new, as is a chapter on advanced data handling. The introductory chapter has been extended and reorganized as two chapters. Exercises have been revised and answers are now provided in an Appendix.

Peter Dalgaard is associate professor at the Department of Biostatistics at the University of Copenhagen and has extensive experience in teaching within the PhD curriculum at the Faculty of Health Sciences. He has been a member of the R Core Team since 1997.


You can download in these links:

link 1

link 2

link 3


To download language R , you shoud visit

http://www.r-project.org/

martes, 7 de diciembre de 2010

Simple OPENGL con CUDA

Aca les dejo un simple ejemplo con OpenGL es el mismo de la pagina de CUDA pero no tiene referenciado nada de sus liberias CUTIL es solo CUDA :). Para compilarlo ejecuta esto "nvcc mesh.cu -o mesh -lglut -lGLEW"

Espero les ayude.





/*
* Copyright 1993-2010 NVIDIA Corporation. All rights reserved.
*
* Please refer to the NVIDIA end user license agreement (EULA) associated
* with this source code for terms and conditions that govern your use of
* this software. Any use, reproduction, disclosure, or distribution of
* this software and related documentation outside the terms of the EULA
* is strictly prohibited.
*
*/

/*
This example demonstrates how to use the Cuda OpenGL bindings to
dynamically modify a vertex buffer using a Cuda kernel.

The steps are:
1. Create an empty vertex buffer object (VBO)
2. Register the VBO with Cuda
3. Map the VBO for writing from Cuda
4. Run Cuda kernel to modify the vertex positions
5. Unmap the VBO
6. Render the results using OpenGL

Host code
*/

// includes, system
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

// includes, GL
#include <GL/glew.h>


#include <GL/glut.h>

// includes
#include <cudaGL.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <vector_types.h>
#include <cuda_gl_interop.h>
#include <unistd.h>




////////////////////////////////////////////////////////////////////////////////
// constants
const unsigned int window_width = 512;
const unsigned int window_height = 512;

const unsigned int width = 256;
const unsigned int height = 256;

// vbo variables
GLuint vbo;
struct cudaGraphicsResource *cuda_vbo_resource;

// mouse controls
int mouse_old_x, mouse_old_y;
int mouse_buttons = 0;
float rotate_x = 0.0, rotate_y = 0.0;
float translate_z = -3.0;

float time1=0.0;


struct cudaGraphicsResource* positionsVBO_CUDA;

/*blocks ands grids*/
dim3 dimBlock(8,8,1);
dim3 dimGrid(width/dimBlock.x,height/dimBlock.y,1);
float4* positions;

////////////////////////////////////////////////////////////////////////////////
//! Keyboard events handler
////////////////////////////////////////////////////////////////////////////////
void keyboard(unsigned char key, int /*x*/, int /*y*/)
{
switch(key) {
case(27) :
exit(0);
break;
}
}

////////////////////////////////////////////////////////////////////////////////
//! Mouse event handlers
////////////////////////////////////////////////////////////////////////////////
void mouse(int button, int state, int x, int y)
{
if (state == GLUT_DOWN) {
mouse_buttons |= 1<<button;
} else if (state == GLUT_UP) {
mouse_buttons = 0;
}

mouse_old_x = x;
mouse_old_y = y;
glutPostRedisplay();
}

void motion(int x, int y)
{
float dx, dy;
dx = x - mouse_old_x;
dy = y - mouse_old_y;

if (mouse_buttons & 1) {
rotate_x += dy * 0.2;
rotate_y += dx * 0.2;
} else if (mouse_buttons & 4) {
translate_z += dy * 0.01;
}

mouse_old_x = x;
mouse_old_y = y;
}


void deleteVBO()
{
cudaGraphicsUnregisterResource(positionsVBO_CUDA);
glDeleteBuffers(1,&vbo);
}

__global__ void createVertices(float4* positions,float time1,unsigned int width,unsigned int height)
{
unsigned int x=blockIdx.x*blockDim.x+threadIdx.x;
unsigned int y=blockIdx.y*blockDim.y+threadIdx.y;

//calculae uv coordinates
float u=x/(float)width;
float v=y/(float)height;

u=u*2.0f-1.0f;
v=v*2.0f-1.0f;

//calculate simple sine wave pattern

float freq=4.0f;
float w=sinf(u*freq+time1)*cosf(v*freq+time1)*0.5f;

//write pòsitions
positions[y*width+x]=make_float4(u,w,v,1.0f);
}


/////////////////////////////////////////////////////////////////////////////////
//Display
////////////////////////////////////////////////////////////////////////////////


void init()
{

cudaGraphicsMapResources(1,&positionsVBO_CUDA,0);
size_t num_bytes;
cudaGraphicsResourceGetMappedPointer((void**)&positions,&num_bytes,positionsVBO_CUDA);




}
void display()
{

//Execute kernel
createVertices<<<dimGrid,dimBlock>>>(positions,time1,width,height);

cudaFree(&positions);
//unmap buffer object

cudaGraphicsUnmapResources(1,&positionsVBO_CUDA,0);

//render from buffer object
glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();

glTranslatef(0.0, 0.0, translate_z);
glRotatef(rotate_x, 1.0, 0.0, 0.0);
glRotatef(rotate_y, 0.0, 1.0, 0.0);

glBindBuffer(GL_ARRAY_BUFFER,vbo);
glVertexPointer(4,GL_FLOAT,0,0);
glEnableClientState(GL_VERTEX_ARRAY);

glColor3f(1.0,0.0,0.0);
glDrawArrays(GL_POINTS,0,width*height);
glDisableClientState(GL_VERTEX_ARRAY);

//swap buffers

glutSwapBuffers();
glutPostRedisplay();


time1+=0.01;
usleep(1000); //sleel 1000 us
}


////////////////////////////////////////////////////////////////////////////////
// Program main
////////////////////////////////////////////////////////////////////////////////

int main(int argc, char** argv)
{

/*paramets to init opengl*/
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
glutInitWindowSize(window_width, window_height);
glutCreateWindow("Cuda GL ");



// initialize necessary OpenGL extensions
glewInit();
if (! glewIsSupported("GL_VERSION_2_0 ")) {
fprintf(stderr, "ERROR: Support for necessary OpenGL extensions missing.");
fflush(stderr);
return -1;
}

// default initialization
glClearColor(0.0, 0.0, 0.0, 1.0);
glDisable(GL_DEPTH_TEST);

// set viewport
glViewport(0, 0, window_width, window_height);

// projection
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
gluPerspective(60.0, (GLfloat)window_width / (GLfloat) window_height, 0.1, 10.0);


//set device NVIDIA

cudaGLSetGLDevice(0);




//create buffer object and register with cuda
glGenBuffers(1,&vbo);
glBindBuffer(GL_ARRAY_BUFFER,vbo); //bind object with vbo
unsigned int size=width*height*4*sizeof(float);
glBufferData(GL_ARRAY_BUFFER,size,0,GL_DYNAMIC_DRAW);
glBindBuffer(GL_ARRAY_BUFFER,0);
cudaGraphicsGLRegisterBuffer(&positionsVBO_CUDA,vbo,cudaGraphicsMapFlagsWriteDiscard);
init(); //map object with cuda

glutDisplayFunc(display);
glutKeyboardFunc(keyboard);
glutMotionFunc(motion);
glutMouseFunc(mouse);

glutMainLoop();
return 0;
}



Unas imagenes de como salio ;)


viernes, 19 de noviembre de 2010

Instalar CUDA GPU NVidia en Ubuntu 10.10

En este tutorial en español describere como pude instalar CUDA Nvidia teniendo una tarjeta de modelo 9500GT. El sistema operativo en el que instale fue Ubuntu 10.10.

Paso 1: (Descargando CUDA)
De este link descarga los siguientes paquetes:
Developer Driver for Linux (260.19.14)
Cuda Toolkit
GPU Computing SDK code sample
Donde descargastes cada uno de los paquetes abre un terminal en esa carpeta y pon : chmod 0777 nombre_paquete

Paso 2: Modifiquemos y removemos nvidia por defecto

En un terminal como modo superusuario hagamos lo sgte:

#nano /etc/modprobe.d/blacklist.conf :
blacklist vga16fb
blacklist rivafb
blacklist nvidiafb
blacklist rivatv

Agregemos las siguientes lineas al final del archivo:
Ejecutemos lo sgte:
#apt-get purge nvidia-*
#reboot -n
Paso 3: Modo X , logueo y parada de servicio gdm
ctrl+alt+f5

nos logueamos y ejecutamos $ sudo service gdm stop

Paso 4: Instalando devdriver
$ sudo sh devdriver_3.2_linux_32_260.19.21.run
Aceptamos todas las licencias
La compatibilidad de OpenGl acepta "Yes"
y preguntara si nvidia-xconfig se desee ejecutar automaticamente para la configuración del archivo X . Responde "Yes"

Paso 5: Instalemos el toolkit de CUDA

Una vez instalado el driver de Nvidia , procedemos a instalar el cudatoolkit. Ejecutando lo siguiente:

$ sudo sh cudatoolkit_3.2.16_linux_32_ubuntu10.04.run
Deja la ruta por defecto , presiona Enter.

Paso 6: Modifiquemos variables del ambiente

Agrega variables del sistema operativo. En un terminal haz lo siguiente:
$sudo nano /etc/environment
Agregemos ":usr/local/cuda/bin" antes de el ultimo simbolo " . Guarda el archivo y salimos del editor

En el archivo .bashrc q esta en tu home. Agregamos lo siguiente:
export CUDA_HOME="/usr/local/cuda"
export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}":${CUDA_HOME}/lib
export PATH=${CUDA_HOME}/bin:${PATH}


# nano /etc/ld.so.conf.d/cuda.conf (creamos nuevo archivo)
/usr/local/cuda/lib64
/usr/local/cuda/lib

Recargamos LD_LIBRARY_PATH

#ldconfig
Paso 7: Instalamos "GPU Computing SDK"

Instalemos gpucomputing : $ sh gpucomputingsdk_3.2.16_linux.run
Enter para instalar el path por defecto.
Dirigimos la ruta de CUDA en este caso esta por defecto en /usr/local/cuda . Presionamos enter para confirmar

Instalamos los paquetes necesarios
$ sudo apt-get install g++ freeglut3-dev libxi-dev
Rompemos el enlace de libGL.so
$ sudo rm -f /usr/lib/libGL.so
$ sudo ln -s /usr/lib/libGL.so.260.19.21 /usr/lib/libGL.so
Creamos un enlaze para libXmu

$sudo ln -s /usr/lib/libXmu.so.6 /usr/lib/libXmu.so

Compilemos GPU COmputing

$ cd ~/NVIDIA_GPU_Computing_SDK/C
$ make

Para encontrar los ejecutable se ubican en ~/NVIDIA_GPU_Computing_SDK/C/bin/linux/release . Uno de sus ejemplos es esta simulacion del humo de cigarro.





Espero les sirva este post!!! xD.




sábado, 13 de noviembre de 2010

Julia SeT in CPU and GPU

Reading the boook cuda by example , find a problem simple "Julia Set" (definition of Julia Set)


CPU



/*
* Copyright 1993-2010 NVIDIA Corporation. All rights reserved.
*
* NVIDIA Corporation and its licensors retain all intellectual property and
* proprietary rights in and to this software and related documentation.
* Any use, reproduction, disclosure, or distribution of this software
* and related documentation without an express license agreement from
* NVIDIA Corporation is strictly prohibited.
*
* Please refer to the applicable NVIDIA end user license agreement (EULA)
* associated with this source code for terms and conditions that govern
* your use of this NVIDIA software.
*
*/


#include "../common/book.h"
#include "../common/cpu_bitmap.h"

#define DIM 1000

struct cuComplex {
float r;
float i;
cuComplex( float a, float b ) : r(a), i(b) {}
float magnitude2( void ) { return r * r + i * i; }
cuComplex operator*(const cuComplex& a) {
return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
}
cuComplex operator+(const cuComplex& a) {
return cuComplex(r+a.r, i+a.i);
}
};

int julia( int x, int y ) {
const float scale = 1.5;
float jx = scale * (float)(DIM/2 - x)/(DIM/2);
float jy = scale * (float)(DIM/2 - y)/(DIM/2);

cuComplex c(-0.8, 0.156);
cuComplex a(jx, jy);

int i = 0;
for (i=0; i<200; i++) {
a = a * a + c;
if (a.magnitude2() > 1000)
return 0;
}

return 1;
}

void kernel( unsigned char *ptr ){
for (int y=0; y<DIM; y++) {
for (int x=0; x<DIM; x++) {
int offset = x + y * DIM;

int juliaValue = julia( x, y );
ptr[offset*4 + 1] = 255 * juliaValue;
ptr[offset*4 + 0] = 0;
ptr[offset*4 + 2] = 0;
ptr[offset*4 + 3] = 255;
}
}
}

int main( void ) {
CPUBitmap bitmap( DIM, DIM );
unsigned char *ptr = bitmap.get_ptr();

kernel( ptr );

bitmap.display_and_exit();
}


GPU





/*
* Copyright 1993-2010 NVIDIA Corporation. All rights reserved.
*
* NVIDIA Corporation and its licensors retain all intellectual property and
* proprietary rights in and to this software and related documentation.
* Any use, reproduction, disclosure, or distribution of this software
* and related documentation without an express license agreement from
* NVIDIA Corporation is strictly prohibited.
*
* Please refer to the applicable NVIDIA end user license agreement (EULA)
* associated with this source code for terms and conditions that govern
* your use of this NVIDIA software.
*
*/


#include "../common/book.h"
#include "../common/cpu_bitmap.h"

#define DIM 1000

struct cuComplex {
float r;
float i;
cuComplex( float a, float b ) : r(a), i(b) {}
__device__ float magnitude2( void ) {
return r * r + i * i;
}
__device__ cuComplex operator*(const cuComplex& a) {
return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
}
__device__ cuComplex operator+(const cuComplex& a) {
return cuComplex(r+a.r, i+a.i);
}
};

__device__ int julia( int x, int y ) {
const float scale = 1.5;
float jx = scale * (float)(DIM/2 - x)/(DIM/2);
float jy = scale * (float)(DIM/2 - y)/(DIM/2);

cuComplex c(-0.8, 0.156);
cuComplex a(jx, jy);

int i = 0;
for (i=0; i<200; i++) {
a = a * a + c;
if (a.magnitude2() > 1000)
return 0;
}

return 1;
}

__global__ void kernel( unsigned char *ptr ) {
// map from blockIdx to pixel position
int x = blockIdx.x;
int y = blockIdx.y;
int offset = x + y * gridDim.x;

// now calculate the value at that position
int juliaValue = julia( x, y );
ptr[offset*4 + 0] = 255 * juliaValue;
ptr[offset*4 + 1] = 0;
ptr[offset*4 + 2] = 0;
ptr[offset*4 + 3] = 255;
}

// globals needed by the update routine
struct DataBlock {
unsigned char *dev_bitmap;
};

int main( void ) {
DataBlock data;
CPUBitmap bitmap( DIM, DIM, &data );
unsigned char *dev_bitmap;

HANDLE_ERROR( cudaMalloc( (void**)&dev_bitmap, bitmap.image_size() ) );
data.dev_bitmap = dev_bitmap;

dim3 grid(DIM,DIM);
kernel<<<grid,1>>>( dev_bitmap );

HANDLE_ERROR( cudaMemcpy( bitmap.get_ptr(), dev_bitmap,
bitmap.image_size(),
cudaMemcpyDeviceToHost ) );

HANDLE_ERROR( cudaFree( dev_bitmap ) );

bitmap.display_and_exit();
}




Really GPU is very interesting , why? The future of computing is parallel computing!!!!!.


See you :P

Get Properties of Device GPU //Obten propiedades de tu dispostivo GPU

Si deseas obtener información de tu tarjeta o tarjetas lo puedes hacer de la siguiente manera:






#include <iostream>

using namespace std;

int main()
{

cudaDeviceProp prop;

int count;
cudaGetDeviceCount(&count); /*count number of devices*/
for(int i=0;i<count;i++)
{
cudaGetDeviceProperties(∝,i); /*get properties of device*/
cout<<"General Information for device----"<<i<<endl;
cout<<"Name "<<prop.name<<endl;
cout<<"Compute capability "<<prop.major<<". "<<prop.minor<<endl;
cout<<"Clock Rate "<<prop.clockRate<<endl;
cout<<"Device copy overlap: ";
if(prop.deviceOverlap)
cout<<"Enable\n";
else cout<<"Disabled\n";

if(prop.kernelExecTimeoutEnabled)
cout<<"Enable\n";
else cout<<"Disabled\n";


cout<<"---Memory Information for device "<<i<<endl;
cout<<"Total global Mem: "<<prop.totalGlobalMem<<endl;
cout<<"Total COnstant Mem: "<<prop.totalConstMem<<endl;
cout<<"Max mem pitch: "<<prop.memPitch<<endl;
cout<<"Texture ALignment "<<prop.textureAlignment;
cout<<"\n -------MP Information for devices "<<i;
cout<<" \nMultiprocessor count : "<<prop.multiProcessorCount<<endl;

cout<<"Shared mem per mp: "<<prop.sharedMemPerBlock<<endl;
cout<<"Registers per mp: "<<prop.regsPerBlock<<endl;
cout<<"Threads in warp: "<<prop.warpSize<<endl;
cout<<"Max Threads per block "<<prop.maxThreadsPerBlock<<endl;
cout<<"Max Threads dimensions: ("<<prop.maxThreadsDim[0]<<","<<prop.maxThreadsDim[1]
<<","<<prop.maxThreadsDim[2]<<")\n";

cout<<"Max grid dimensions: ("<<prop.maxGridSize[0]<<","<<prop.maxGridSize[1]
<<","<<prop.maxGridSize[2]<<")\n";
}
return 0;
}


see you!!! xD

Multiplicación de Matricez Multibloque con CUDA

Ahora multiplicación de Matrizes , pero usando un nuevo concepto TILES , en el cual podemos hacer multiplicación de matrices multibloque.









Acontinuación el sgte codigo en el que aplica ese concepto:




#include <stdio.h>
#include <assert.h>
#include <iostream>

#define TILE_WIDTH 4
using namespace std;

// Simple utility function to check for CUDA runtime errors
void checkCUDAError(const char *msg);

// Part 3 of 5: implement the kernel
__global__ void MatrixTiles(int *d_a ,int *d_b,int *d_c,int width )
{
int row=blockIdx.x*TILE_WIDTH+threadIdx.x;
int col=blockIdx.y*TILE_WIDTH+threadIdx.y;


int pvalue=0;
//for each computes one element of the block sub-matrix
for(int i=0;i<width;i++)
{
pvalue+=(d_a[row*width+i]*d_b[i*width+col]);
}
d_c[row*width+col]=pvalue;
}

////////////////////////////////////////////////////////////////////////////////
// Program main
////////////////////////////////////////////////////////////////////////////////
int main( int argc, char** argv)
{
// pointer for host memory
int *h_a,*h_b,*h_c;

// pointer for device memory
int *d_a,*d_b,*d_c;

// define grid and block size
int numBlocks = 4;
int numThreadsPerBlock = 16;

// Part 1 of 5: allocate host and device memory
size_t memSize = numBlocks * numThreadsPerBlock * sizeof(int);
h_a = (int *) malloc(memSize);
cudaMalloc((void**) &d_a,memSize );


h_b = (int *) malloc(memSize);
cudaMalloc((void**) &d_b,memSize );

h_c = (int *) malloc(memSize);
cudaMalloc((void**) &d_c,memSize );



for(int n=0;n<numThreadsPerBlock*numBlocks;n++)
{
h_a[n]=h_b[n]=1;
}

// Part 2 of 5: configure and launch kernel
dim3 dimGrid( numBlocks/2,numBlocks/2 );
dim3 dimBlock( numThreadsPerBlock/TILE_WIDTH,numThreadsPerBlock/TILE_WIDTH);




// check if kernel execution generated an error
checkCUDAError("kernel execution");

// Part 4 of 5: device to host copy
cudaMemcpy( d_a,h_a ,memSize,cudaMemcpyHostToDevice);
cudaMemcpy( d_b,h_b ,memSize,cudaMemcpyHostToDevice);


MatrixTiles<<< dimGrid ,dimBlock >>>(d_a,d_b,d_c,TILE_WIDTH*2);
// block until the device has completed
cudaThreadSynchronize();
cudaMemcpy( h_c,d_c ,memSize,cudaMemcpyDeviceToHost);

// Check for any CUDA errors
checkCUDAError("cudaMemcpy");

for(int i=0;i<8;i++)
{
for(int j=0;j<8;j++)
cout<<h_c[8*i+j]<<" ";
cout<<endl;
}

// free device memory
cudaFree(d_a);

// free host memory
free(h_a);

// If the program makes it this far, then the results are correct and
// there are no run-time errors. Good work!
printf("Correct!\n");

return 0;
}

void checkCUDAError(const char *msg)
{
cudaError_t err = cudaGetLastError();
if( cudaSuccess != err)
{
fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString( err) );
exit(-1);
}
}




El sgte post será mejorando esto. Los Tiles se usan porque cada bloque en CUDA tienen un tamaño que es de 512 Threads.

lunes, 8 de noviembre de 2010

Mi primer programa en GPU -CUDA


Despues de unas leidas para unos trabajos que hare mi primer programa en GPU aún es una simple multiplicación de matrizes con Hilos en CUDA. Las versiones mejoradas ya las subire.




#include <stdio.h>
#include <assert.h>
#include <iostream>


using namespace std;

// Simple utility function to check for CUDA runtime errors
void checkCUDAError(const char *msg);

// Part 3 of 5: implement the kernel
__global__ void MatrixSimple(int *d_a ,int *d_b,int *d_c,int width )
{
int row=blockIdx.x*width+threadIdx.x;
int col=blockIdx.y*width+threadIdx.y;
int pvalue=0;
//for each computes one element of the block sub-matrix
for(int i=0;i<width;i++)
{
pvalue+=(d_a[row*width+i]*d_b[i*width+col]);
}
d_c[row*width+col]=pvalue;
}

////////////////////////////////////////////////////////////////////////////////
// Program main
////////////////////////////////////////////////////////////////////////////////
int main( int argc, char** argv)
{
// pointer for host memory
int *h_a,*h_b,*h_c;

// pointer for device memory
int *d_a,*d_b,*d_c;

// define grid and block size
int numBlocks = 1;
int numThreadsPerBlock = 64;

// Part 1 of 5: allocate host and device memory
size_t memSize = numBlocks * numThreadsPerBlock * sizeof(int);
h_a = (int *) malloc(memSize);
cudaMalloc((void**) &d_a,memSize );
h_b = (int *) malloc(memSize);
cudaMalloc((void**) &d_b,memSize );
h_c = (int *) malloc(memSize);
cudaMalloc((void**) &d_c,memSize );
for(int n=0;n<numThreadsPerBlock*numBlocks;n++)
{
h_a[n]=h_b[n]=1;
}

// Part 2 of 5: configure and launch kernel
dim3 dimGrid( numBlocks );
dim3 dimBlock( numThreadsPerBlock/8,numThreadsPerBlock/8);


// check if kernel execution generated an error
checkCUDAError("kernel execution");

// Part 4 of 5: device to host copy
cudaMemcpy( d_a,h_a ,memSize,cudaMemcpyHostToDevice);
cudaMemcpy( d_b,h_b ,memSize,cudaMemcpyHostToDevice);
MatrixSimple<<< dimGrid ,dimBlock >>>(d_a,d_b,d_c,numThreadsPerBlock/8);
// block until the device has completed
cudaThreadSynchronize();
cudaMemcpy( h_c,d_c ,memSize,cudaMemcpyDeviceToHost);

// Check for any CUDA errors
checkCUDAError("cudaMemcpy");

for(int i=0;i<8;i++)
{
for(int j=0;j<8;j++)
cout<<h_c[8*i+j]<<" ";
cout<<endl;
}

// free device memory
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);

// free host memory
free(h_a);
free(h_b);
free(h_c);

// If the program makes it this far, then the results are correct and
// there are no run-time errors. Good work!
printf("Correct!\n");

return 0;
}

void checkCUDAError(const char *msg)
{
cudaError_t err = cudaGetLastError();
if( cudaSuccess != err)
{
fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString( err) );
exit(-1);
}
}


Hasta la prox!!!

domingo, 24 de octubre de 2010

Nuestro Primer ACM-ICPC en Trujillo-Perú part2 (UNT)

parte 2: El inicio del ACM-ICPC en la UNT


Comenzar algo siempre es lo complicado en toda la vida por ejemplo: ingresar a la U, lanzarte a tu primera flaka, y creo que esta lista no terminaría.

Como en mis post anterior escribí que hace como 2 años la UNT quería ir a este concurso , pero esa gente ya salio y solo quedamos 2. Después de todo el tema organizacional y al haber obtenido la sede del icpc en mi UNT (mi U ) era el gran problema tener gente dispuesta a participar ese fue un problema casi NP xD.


Ya a comienzos de año entre salón por salón a invitar gente para que participen pero nada a las reunión solo fueron :S 4 gatos y poco a poco solo quedamos 2 como al comienzo :( , de allí seguimos como 4 a 6 semanas así , en 2 semanas mas replanteamos todo y buscamos gente mas directa y de confianza y con habilidades :D aun comienzo se juntaron como 28 waooo con esta gente si la hacemos dije :\ .

Pero creo que ellos al no encontrar coach ni nada por el estilo no vieron un trabajo serio a las finales solo resultamos 7 :| , alli entrenábamos solo sin nadie que nos ayude ("que novedad en la UNT") , hasta llego un tiempo en el que nadie entrenaba :S una tremenda preocupación esto porque ponía en riesgo el acm-icpc :S en Trujillo.

Termino ese semestre y no fueron los entrenamientos como hubiéramos querido :( , así que se me ocurrió otra loca o brillante idea (por lo general son ideas locas xD) , hacer de asistente de cátedra para este semestre que ya paso hace una semana, con esto podría sacar mas gente pero de ciclos inferiores y que ellos continúen con este legado :P .

Así que me metí a ser asistente aun en mi septimo ciclo ; el curso del que se trataba era sobre análisis y diseño de algoritmos (exacto lo que viene en el acm-icpc) por eso me anime mucho, la idea es buena pero había un problema ahora :( tenia que reducirme cursos para poder hacer esto , no se si tome una buena o mala decisión pero deje 1 curso y solo lleve 5 cursos , para que esto corra .

Asi me metí a hacer asistencia en dicho curso y bueno , pobres mis alumnos :S los traume las primeras clases, problemas en ingles y encima los jueces que los rechazaba xD . Pobres ya lloraban, tanto fue que hasta emails jodiendo me mandaron ¬¬!! cosa que no me agrado , poco tiempo después S: tuve problemas con creditaje y pues me meti a 2 cursos mas en total tenia 7 cursos+asistente= 8 cursos . Waooo mi vida se volvió extremadamente horrible xD . Y eso sin contar que estuve en la organización del congreso de la SPC

Ya de alli pude seleccionar a 2 buenos alumnos y para que conozcan las ligas mayores de programación :P .

Los otros 2 miembros de mi team (Leandro y Karla) también fueron asistentes de profesores xD , uno en algoritmos (1 ciclo) y otro programación lógica.

Lo mejor de esto es que mi amigo Leando les mostró también los jueces para entrenar y ya tenemos cachimbos entrenando :') csm!!! no esperaba esto ni lo soñaba aun xD, esto es excelente para futuras participaciones.

De toda la gente que nos acompaño en el ACM-ICPC realmente les estoy como 1/0 (indeterminadamente) agradecido por su participación en este concurso , por darle un tiempo a esto, yo se que a pesar de muchas variables que jugaron en contra nuestra incluso siendo todos organizadores del evento xD, y haber tenido una semana antes un congreso .

A pesar de todo lo malo que nos paso y sin tener el apoyo directo de un coach ni nadie que nos ayudara en esto supimos manejarlo , hubiéramos querido hacer mas pero para ser un primer año los resultados no son malos .


No dare mas detalles de todo lo que paso en la conformación de nuestro equipos xD pero de que si entrenamos si lo hicimos , solos pero alli estuvimos ;).


Ahora si hemos dado el primer paso , esto fue lo mas complicado de todo comenzar!!! , alli en toda nuestra terquedad y ganas de que nuestra escuela participe en un evento de nivel internacional como este :'). Y gracias a nuestra querida UNT que nos financio para este evento al menos en esto nos ayudaron :D .


Ahora se viene todo un año duro de entrenamiento :D y demás cosas ;) , espero antes de morirme ver a un equipo de mi Universidad Nacional de Trujillo en una final mundial del ACM-ICPC (y si es en el podio de ganador mejor ;)), así como fui y fuimos de tercos , lo lograremos algún día gogogogogo!!!

La frase de todo ese tiempo: "No esperes que tu Universidad te de; tu dale a tu Universidad".



PD: Alli una foto con toda la gente del team UNT y con el profe Corrado que ayudo enormemente ;).



Nuestro Primer ACM-ICPC en Trujillo-Perú part1


Parte 1: La Organización

Como ya todos sabemos este año se realizo el primer acm-icpc en Trujillo , algo nunca antes visto en nuestra ciudad , ni siquiera pensado hace unos años que esto se daría. Todo comenzó en diciembre del 2009, yo ya sabia que era el acm-icpc y hace como 2 años atras que queríamos participar en este concurso , pero los costos de viaje + hospedaje+ inscripción , difícilmente hubiéramos podido pagar porque somos de una Universidad Nacional.

Como tuvimos referencia que la UPAO había ido el 2009 a Lima , y también sabia que la UPN ya había participado , se me ocurrió la brillante o alocada idea de traer el concurso a Trujillo. Comenze haciendo el plan de hormiga , visite al director de escuela de sistemas de UPAO y le comente si el nos apoyaría para traer esta sede , me comento que el ya había conversado con Viktor Knelokovf(Representante de acm-icpc en Perú) para traerlo dicho evento , pero no obtuvo buenas respuestas , coordinamos bien y el extendió todo su apoyo a esto.

Acto seguido visitar la UPN para coordinar si nos apoyaría en este camino xD, también me comentaron que trataron de traer dicho evento solos pero nada de nada , y también me extendió su apoyo a dicho evento.

Y como última parte en este Enero converse con el director de mi escuela y también acepto (como siempre xD). Emails por aca por aya , converse con el Dr. Cesa Beltrán si nos podía ayudar con esto para seguir su esquema del concurso como ellos lo llevaron dicho concurso a Arequipa.

Ya para febrero primeros días tuvimos una reunión UPN,UPAO y UNT ; cosa rara que se junten para algo, allí conversamos sobre su apoyo a dicho evento y la UNT se eligió como sede :D . Después de emails con Viktor vimos que implicaba otros gastos y con 3 universidades no cubríamos gasto :\ . Asi que llamamos a la UCV ¬¬ , acepto participar y todo un mes despues de eso otra reunión con todo ellos para asegurar dicha participación.

Ya unos pocos meses despues Viktor nos acepto ser sede :D , bueno le enviamos una carta de compromiso en el que todas las U's se comprometían a participar para asegurar dicho evento :D .

Ya la organización local todos fuimos contestant xD pero weno no querían apoyar profes asi que ya pes asumimos el reto.


Ya no dare mas detalle , pero este evento en verdad es extremandamente complicado hay muchos problemas siempre :| toda esta semana que paso andaba al borde del suicidio de tanto stress porque todo tenia que salir bien, y menos mal que salio todo ook. Por alli que a comienzo de semana la UCV declino su participación y nos ajusto mucho todos nuestros planes. Y cosas asi paso toda está semana.

Y al siguiente año , como esta sede rotará para que el mundo tenga un equilibrio , quedaría entre la UPAO o UPN o quiza otra U que asuma este reto.

Nos reencontramos el siguiente año y para esa vez la UNT no cometera la burrada de que sus contestant sean tambien organizadores.

Y agradecer a Viktor por la sede :D , Alfredo Paz por el gran apoyo y al prof. Corrado que nos ayudo enormemente en este concurso , y a todos los demás de la UNT (prof vilma,...etc)

Adios y mi frase de toda esta semana "Lo mejor que me puede pasar es lo peor".

PD: Abajo nuestro sipan programador xD





sábado, 27 de marzo de 2010

Google File System (GFS)

Google quién no lo usa y se da cuenta que es muy rápido, es que google como una verdadera empresa de Ciencia de la Computación en realidad sabe innovar como es en cada aplicación que desarrollada , es verdaderamente fascinante lo que gracias a google avanzo la computación que hoy en dia es el gran competidor de Microsoft =).

Bueno la historia anterior casi la mayoria lo sabe. Lo que les comentaré es sobre el Google File System ; es tanto lo que avanzo Google que puede manejar petabytes de almacenamiento algo increible en verdad q hace unos 6 a 8 años atras ni pensabamos que llegariamos a esa cantidad de datos.

Su modelo el GFS es tan interesante que practicamente tiro al tacho todas las teorias de SFD(System File Distributed) pues no tiene comparación ni con NFS ("Networking File System") ni con AFS ("Andres File System") .

En la última revista de ACM , tiene un articulo interesante que nos comenta sobre como evoluciono su sistema de archivos en si es una conversación en la que los dos autores empiezan a narrar, lo anecdótico es que iban prácticamente chancando cada vez que agregaban un nuevo app de google.

En este link se encuentra la entrevista y en este otro en version pdf.

Nos vemos :D

sábado, 20 de marzo de 2010

Configurar OpenMesh en Linux / Configure Open Mesh in Linux



Este semestre que llevo computación gráfica 2, la verdad deja mucho que desear, mis ideas pa este ciclo era hacer mallas, tanto fue mi obsesion que lei CGAL un monton y nada no me salio :S.
Asi que ya de tanto buscar encontre OpenMesh =), una libreria para manipulación de mallas y polygonos.

Ahora describere como lo instale en mi computador :
  • Descargas la librería desde aquí .
  • Descomprimelo con el gestor de archivadores en el lugar que desees.
  • Ahora ejecuta cmake-gui como superusuario. Y rutea desde donde descomprimistes y a donde se compilara. Presiona configure y esperas q termine con exito y cierras.


  • En la carpeta donde se compilo. Ingresa en la carpeta buil. Entonces abre un terminal como superusuario y ponerle "cmake .."


  • Ahora regresa a su raiz y compila com el comando make.


  • Esperamos un momento que termine.


  • Terminado nos dirigimos a la carpeta /bin donde se compilo y ejecutamos QtViewer , para probarlo descargamos el dragon de standford de este link . Abirmos el archivo y esperamos que carge esto demorara un poco :D.



  • Espero seguir con esto de computación gráfica aunque tenga un pesimo profesor. Hasta otro post

    viernes, 12 de marzo de 2010

    Code jam Minimum Scalar Product/Menor producto escalar



    //Problema: Menor producto escalar
    //Desarrollado por: Marks Calderon Niquin

    #include <iostream>
    #include <vector>
    #include <fstream>
    #include <algorithm>
    #include <cstdio>
    #include <cstdlib>


    using namespace std;



    int main()

    {
    ifstream in;
    ofstream out;
    in.open("A-large-practice.in",ifstream::in);
    out.open("salida.out");
    int n,n1;
    vector<long long> p,q;
    long long temp;
    in>>n;
    for(int i=0;i<n;i++)
    {
    in>>n1;

    for(int k=0;k<n1;k++)
    {
    in>>temp;
    p.push_back(temp);
    }
    for(int k=0;k<n1;k++)
    {
    in>>temp;
    q.push_back(temp);
    }
    sort(p.begin(),p.end()); //sort in two vectors
    sort(q.begin(),q.end());
    temp=0;
    for(int l=0;l<n1;l++) //multiplicate
    temp+=(p[l]*q[q.size()-1-l]);
    out<<"Case #"<<(i+1)<<": "<<temp<<endl;
    p.clear();
    q.clear();
    }
    in.close();
    out.close();


    return 0;

    }

    100-3n+1

    The famous problem 3n+1 is a problem NP. Try to improve my time


    /*

    Judge : Uva
    Number: 100
    Name : 3n+1
    Time: 0.028
    */

    #include <stdio.h>
    #include <stdlib.h>
    #define MAX 1000001
    #define SWAP(a, b) { a ^= b;b ^= a;a ^= b;} /*swapping*/

    unsigned long conta[MAX];

    unsigned long tres(unsigned long n)
    {
    if(n<MAX)
    {
    if(conta[n]==0) conta[n]=(n&1)?1+tres(3*n+1):1+tres(n>>1);
    return conta[n];
    }
    else
    {
    if(n%2) return 2+ tres((3*n+1)>>1); /*how n is > MAX */
    else return 1+ tres(n>>1);
    }

    }
    int main()
    {

    unsigned long a,b,c,d,maximo,i,temp;
    for(i=0;i<=MAX;i++) //init table
    conta[i]=0;
    conta[1]=1;
    while(scanf("%lu%lu",&a,&b)==2)
    {
    c=a;
    d=b;
    if(a>b) SWAP(a,b);
    maximo=0;
    for(i=a;i<=b;i++)
    {
    if(conta[i]==0)
    conta[i]=tres(i);
    if(maximo<conta[i])
    maximo=conta[i];
    }
    printf("%lu %lu %lu\n",c,d,maximo);
    }
    return 0;
    }