Re: [Help-glpk] Python glpk

[Carlos Herrera]

El mar, 05-06-2007 a las 09:33 -0700, J.-O. Moussafir escribió:
> Hello,
> 
> I've come across segfault using Python-glpk. Does anyone have had some
experience with it ?
> 
> Jacques-Olivier
> 
--------------------------------------------------------------------------------------------------------------------------

> I have experience with PULP (A linear programming modeler in Python -
http://www.jeannot.org/~js/code/index.en.html#PuLP).
> This is a program developed in python to use GLPK in python. I send
you an example of "data envelopment analysis" 
> developed using Pulp. 

###############################################################
#!/usr/bin/env python
# @(#) $Jeannot: test1.py,v 1.11 2005/01/06 21:22:39 js Exp $
# Importar PuLP: A python linear programming modeler
from pulp import *
# Importar Scipy: Scientific python
from scipy import *

from string import *
# Importar otras librerias
import Gnuplot, Gnuplot.funcutils,os,time,datetime

# Funciones a utilizar
# Funcion para leer datos
def LeeDat(archivo): # lee datos del archivo
  f = open(archivo, 'r')
  salida = io.read_array(f)
  f.close()
  return salida

# Funcion para generar archivos
def Guarda(variable,nombre):  # Guarda variable en archivo
  archivo = file(nombre, 'w') # Creando archivo
  io.write_array(archivo,variable) # Escribiendo archivo
  archivo.close() # Cerrando archivo

# Comienza el programa
intro = "\n\
                              Data Envelopment Analysis\n\
                                  Version Beta-0.1\n\
                            Carlos Herrera - Andres Vejar\n\
                        DEPARTAMENTO DE INGENIERIA INDUSTRIAL\n\
                           UNIVERSIDAD DE SANTIAGO DE CHILE\n\
                           ::::::::::::::::::::::::::::::::"
print intro
print ""
print "Parametros del modelo"
print "---------------------"
i1 = 19
print "El numero de entidades a evaluar es",i1
print ""
i2 = int(raw_input("Ingrese el numero de meses: "))

# Parametros del modelo
dmu =  i1 # numero de entidades
meses = i2 # numero de meses

# Matriz de ponderaciones de eficiencia
solu0 = zeros((dmu,meses))
solu1 = mat(zeros((dmu,meses)))
solu2 = mat(zeros((dmu,meses)))
for mes in range(meses):
  entrada0 = 'entrada-' + str(mes+1) +'.datos'
  salida0  = 'salida-' + str(mes+1) +'.datos'
  entrada1 = 'entrada-' + str(mes+2) +'.datos'
  salida1  = 'salida-' + str(mes+2) +'.datos'
  datIn0  = LeeDat(entrada0)
  datOut0 = LeeDat(salida0)
  datIn1  = LeeDat(entrada1)
  datOut1 = LeeDat(salida1)
  entr = shape(datIn0)[1]
  sal = shape(datOut0)[1]
  
  
# ANALISIS ENVOLVENTE DE DATOS (DATA ENVELOPMENT ANALYSIS)

# eficiencia T (x,y)^T
  for dato in range(dmu):
    # Un nuevo programa LP objeto prob nombre "hola"
    dea = LpProblem("dea", LpMinimize)
    # Definicion de la variable eficiencia
    eficiencia = LpVariable("eficiencia", 0)
    # Definicion de la variable lambda
    lambd = LpVariable.matrix("lambda", range(dmu), 0)

    # Definicion de funcion Objetivo
    dea += eficiencia

    for item in range(entr):
      dea += eficiencia*datIn0[dato,item] -
lpDot(list(datIn0[:,item]),lambd) >= 0
   
    for item in range(sal):
      dea += -datOut0[dato,item] + lpDot(list(datOut0[:,item]),lambd) >=
0
  
    dea += lpSum(lambd) == 1
    dea.solve()
  
    solu0[dato,mes] = value(dea.objective)

# eficiencia T (x,y)^T+1
  for dato in range(dmu):
    datIn1[dato,:] = datIn0[dato,:]
    datOut1[dato,:] = datOut0[dato,:]
    # Un nuevo programa LP objeto prob nombre "hola"
    dea = LpProblem("dea", LpMinimize)
    # Definicion de la variable eficiencia
    eficiencia = LpVariable("eficiencia", 0)
    # Definicion de la variable lambda
    lambd = LpVariable.matrix("lambda", range(dmu), 0)

    # Definicion de funcion Objetivo
    dea += eficiencia

    for item in range(entr):
      dea += eficiencia*datIn0[dato,item] -
lpDot(list(datIn1[:,item]),lambd) >= 0
   
    for item in range(sal):
      dea += -datOut0[dato,item] + lpDot(list(datOut1[:,item]),lambd) >=
0
  
    dea += lpSum(lambd) == 1
    dea.solve()
  
    solu1[dato,mes] = value(dea.objective)
    
# eficiencia T+1 (x,y)^T
  for dato in range(dmu):
    datIn0[dato,:] = datIn1[dato,:]
    datOut0[dato,:] = datOut1[dato,:]
    # Un nuevo programa LP objeto prob nombre "hola"
    dea = LpProblem("dea", LpMinimize)
    # Definicion de la variable eficiencia
    eficiencia = LpVariable("eficiencia", 0)
    # Definicion de la variable lambda
    lambd = LpVariable.matrix("lambda", range(dmu), 0)

    # Definicion de funcion Objetivo
    dea += eficiencia

    for item in range(entr):
      dea += eficiencia*datIn1[dato,item] -
lpDot(list(datIn0[:,item]),lambd) >= 0
   
    for item in range(sal):
      dea += -datOut1[dato,item] + lpDot(list(datOut0[:,item]),lambd) >=
0
  
    dea += lpSum(lambd) == 1
    dea.solve()
  
    solu2[dato,mes] = value(dea.objective)
    print 'VER ESTO !!!!!', value(dea.objective)
#print 'NORMAL'
#print solu0
#print 'E (T) - (x,y)^T+1'
#print solu1
Guarda(solu2,'resultados1.txt')

# INDICE DE MALMQUIST

malmquist = mat(zeros((dmu,meses)))
EC = mat(zeros((dmu,meses)))
PT = mat(zeros((dmu,meses)))
for mes in range(meses-1):
  for d in range(dmu):
    EC[d,mes+1] = solu0[d,mes+1]/solu0[d,mes]
    PT[d,mes+1] = ((solu1[d,mes]/solu0[d,mes
+1])*(solu0[d,mes]/solu2[d,mes]))**(1.0/2.0)
    malmquist[d,mes+1] = EC[d,mes+1] * PT[d,mes+1]
 
#########################################################################################

Hopefully that serves to you, 

greetings, 

Carlos
> _______________________________________________
> Help-glpk mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-glpk