線形計画法

# -*- coding: UTF-8 -*-
from math import *
import numpy as np

############################################
# クラス LP(線形計画法)
#      coded by Y.Suganuma
############################################

class LP :

	################################
	# コンストラクタ
	#      n1 : 制約条件の数
	#      m1 : 変数の数
	#      z1 : 目的関数の係数
	#      eq_l : 制約条件の左辺
	#      eq_r : 制約条件の右辺
	#      cp1 : 比較演算子
	################################

	def __init__(self, n1, m1, z1, eq_l, eq_r, cp1) :
				# 初期設定
		self.eps = 1.0e-10   # 許容誤差
		self.err = 0   # エラーコード (0:正常終了, 1:解無し)
		self.n   = n1   # 制約条件の数
		self.m   = m1   # 変数の数
		self.a   = np.zeros(self.n)   # 人為変数があるか否か
		self.cp  = cp1   # 比較演算子(-1:左辺<右辺, 0:左辺=右辺, 1:左辺>右辺)
		self.z   = z1   # 目的関数の係数
				# スラック変数と人為変数の数を数える
		self.m_s = 0   # スラック変数の数
		self.m_a = 0   # 人為変数の数
		for i1 in range(0, self.n) :
			if self.cp[i1] == 0 :
				self.m_a += 1
				if eq_r[i1] < 0.0 :
					eq_r[i1] = -eq_r[i1]
					for i2 in rang(0, self.m) :
						eq_l[i1][i2] = -eq_l[i1][i2]
			else :
				self.m_s += 1
				if eq_r[i1] < 0.0 :
					self.cp[i1] = -self.cp[i1]
					eq_r[i1]    = -eq_r[i1]
					for i2 in rang(0, self.m) :
						eq_l[i1][i2] = -eq_l[i1][i2]
				if self.cp[i1] > 0 :
					self.m_a += 1
				# 単体表の作成
					# 初期設定
		self.mm  = self.m + self.m_s + self.m_a
		self.row = np.empty(self.n, np.int)   # 各行の基底変数の番号
		self.s   = np.empty((self.n+1, self.mm+1), np.float)   # 単体表
		for i1 in range(0, self.n+1) :
			if i1 < self.n :
				self.s[i1][0] = eq_r[i1]
				for i2 in range(0, self.m) :
					self.s[i1][i2+1] = eq_l[i1][i2]
				for i2 in range(self.m+1, self.mm+1) :
					self.s[i1][i2] = 0.0
			else :
				for i2 in range(0, self.mm+1) :
					self.s[i1][i2] = 0.0
					# スラック変数
		k = self.m + 1
		for i1 in range(0, self.n) :
			if self.cp[i1] != 0 :
				if self.cp[i1] < 0 :
					self.s[i1][k] = 1.0
					self.row[i1]  = k - 1
				else :
					self.s[i1][k] = -1.0
				k += 1
					# 人為変数
		for i1 in range(0, self.n) :
			if self.cp[i1] >= 0 :
				self.s[i1][k] = 1.0
				self.row[i1]  = k - 1
				self.a[i1]    = 1
				k += 1
					# 目的関数
		if self.m_a == 0 :
			for i1 in range(0, self.m) :
				self.s[self.n][i1+1] = -self.z[i1]
		else :
			for i1 in range(0 , self.m+self.m_s+1) :
				for i2 in range(0 , self.n) :
					if self.a[i2] > 0 :
						self.s[self.n][i1] -= self.s[i2][i1]

	################################
	# 最適化
	#      sw : =0 : 中間結果無し
	#           =1 : 中間結果あり
	#      return : =0 : 正常終了
	#             : =1 : 解無し
	################################

	def optimize(self, sw) :
				# フェーズ1
		if sw > 0 :
			if self.m_a > 0 :
				print("\nphase 1")
			else :
				print("\nphase 2")
		self.opt_run(sw)
				# フェーズ2
		if self.err == 0 and self.m_a > 0 :
					# 目的関数の変更
			self.mm -= self.m_a
			for i1 in range(0, self.mm+1) :
				self.s[self.n][i1] = 0.0
			for i1 in range(0 , self.n) :
				k = self.row[i1]
				if k < self.m :
					self.s[self.n][0] += self.z[k] * self.s[i1][0]
			for i1 in range(0, self.mm) :
				for i2 in range(0, self.n) :
					k = self.row[i2]
					if k < self.m :
						self.s[self.n][i1+1] += self.z[k] * self.s[i2][i1+1]
				if i1 < self.m :
					self.s[self.n][i1+1] -= self.z[i1]
					# 最適化
			if sw > 0 :
				print("\nphase 2")
			self.opt_run(sw)

		return self.err

	################################
	# 最適化(単体表の変形)
	#      sw : =0 : 中間結果無し
	#           =1 : 中間結果あり
	################################

	def opt_run(self, sw) :

		self.err = -1
		while self.err < 0 :
				# 中間結果
			if sw > 0 :
				print()
				self.output()
				# 列の選択(巡回を防ぐため必ずしも最小値を選択しない,Bland の規則)
			q = -1
			for i1 in range(1, self.mm+1) :
				if self.s[self.n][i1] < -self.eps :
					q = i1 - 1
					break
				# 終了(最適解)
			if q < 0 :
				self.err = 0
				# 行の選択( Bland の規則を採用)
			else :
				p   = -1
				k   = -1
				min = 0.0
				for i1 in range(0, self.n) :
					if self.s[i1][q+1] > self.eps :
						x = self.s[i1][0] / self.s[i1][q+1]
						if (p < 0) or (x < min) or ((x == min) and (self.row[i1] < k)) :
							min = x
							p   = i1
							k   = self.row[i1]
					# 解無し
				if p < 0 :
					self.err = 1
					# 変形
				else :
					x           = self.s[p][q+1]
					self.row[p] = q
					for i1 in range(0, self.mm+1) :
						self.s[p][i1] /= x
					for i1 in range(0, self.n+1) :
						if i1 != p :
							x = self.s[i1][q+1]
							for i2 in range(0, self.mm+1) :
								self.s[i1][i2] -= (x * self.s[p][i2])

	################################
	# 単体表の出力
	################################

	def output(self) :

		for i1 in range(0, self.n+1) :
			if i1 < self.n :
				print("x" + str(self.row[i1]+1), end="")
			else :
				print(" z", end="")
			for i2 in range(0, self.mm+1) :
				print(" " + str(self.s[i1][i2]), end="")
			print()

	################################
	# 結果の出力
	################################

	def result(self) :

		if self.err > 0 :
			print("\n解が存在しません")
		else :
			print("\n(", end="")
			for i1 in range(0, self.m) :
				x = 0.0
				for i2 in range(0, self.n) :
					if self.row[i2] == i1 :
						x = self.s[i2][0]
						break
				if i1 == 0 :
					print(x, end="")
				else :
					print(", " + str(x), end="")
			print(") のとき,最大値 " + str(self.s[self.n][0]))

----------------------------------

# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
from function import LP

############################################
# 線形計画法
#      coded by Y.Suganuma
############################################

			# 入力
				# 変数の数と式の数
line = sys.stdin.readline()
x    = line.split()
m    = int(x[0])
n    = int(x[1])
cp   = np.empty(n, np.int)
z    = np.empty(m, np.float)
eq_r = np.empty(n, np.float)
eq_l = np.empty((n, m), np.float)
				# 目的関数の係数
line = sys.stdin.readline()
x    = line.split()
for i1 in range(0, m) :
	z[i1] = float(x[i1])
				# 制約条件
for i1 in range(0, n) :
	line = sys.stdin.readline()
	x    = line.split()
	for i2 in range(0, m) :
		eq_l[i1][i2] = float(x[i2])
	if x[m] == '<' :
		cp[i1] = -1
	elif x[m] == '>' :
		cp[i1] = 1
	else :
		cp[i1] = 0
	eq_r[i1] = float(x[m+1])
			# 実行
lp = LP(n, m, z, eq_l, eq_r, cp)
lp.optimize(1)
			# 結果の出力
lp.result()