#781385 (PyPy3) No.453 製薬会社

提出ソース
結果

問題	No.453 製薬会社
ユーザー	toyuzuko
提出日時	2022-07-25 00:37:53
言語	PyPy3 (7.3.15)
結果	AC
実行時間	79 ms / 2,000 ms
コード長	13,971 bytes
コンパイル時間	618 ms
コンパイル使用メモリ	87,268 KB
実行使用メモリ	71,724 KB
最終ジャッジ日時	2023-09-21 01:41:11
合計ジャッジ時間	2,215 ms
ジャッジサーバーID （参考情報）	judge12 / judge14
このコードへのチャレンジ
（要ログイン）
テストケース

テストケース表示
入力	結果	実行時間実行使用メモリ
testcase_00	AC	77 ms 71,416 KB
testcase_01	AC	78 ms 71,624 KB
testcase_02	AC	76 ms 71,480 KB
testcase_03	AC	77 ms 71,596 KB
testcase_04	AC	76 ms 71,432 KB
testcase_05	AC	78 ms 71,568 KB
testcase_06	AC	78 ms 71,624 KB
testcase_07	AC	77 ms 71,600 KB
testcase_08	AC	79 ms 71,492 KB
testcase_09	AC	78 ms 71,724 KB
testcase_10	AC	78 ms 71,200 KB
testcase_11	AC	77 ms 71,172 KB
testcase_12	AC	77 ms 71,524 KB
権限があれば一括ダウンロードができます
ソースコード

raw source code
EPS = 1e-10

def is_same(x, y):
    if abs(x - y) < EPS:
        return True
    return False

def is_zero(x):
    if abs(x) < EPS:
        return True
    return False

class Matrix():
    def __init__(self, n, m, mat=None):
        self.n = n
        self.m = m
        self.mat = [[0.] * self.m for _ in range(self.n)]
        if mat:
            assert len(mat) == n and len(mat[0]) == m
            for i in range(self.n):
                self.mat[i] = mat[i].copy()

    def is_square(self):
        return self.n == self.m

    def __getitem__(self, key):
        if isinstance(key, slice):
            return self.mat[key]
        else:
            assert key >= 0
            return self.mat[key]

    def id(n):
        res = Matrix(n, n)
        for i in range(n):
            res[i][i] = 1.
        return res

    def __len__(self):
        return len(self.mat)

    def __str__(self):
        return '\n'.join(' '.join(map(str, self[i])) for i in range(self.n))

    def times(self, k):
        res = [[0.] * self.m for _ in range(self.n)]
        for i in range(self.n):
            res_i, self_i = res[i], self[i]
            for j in range(self.m):
                res_i[j] = k * self_i[j]
        return Matrix(self.n, self.m, res)

    def __pos__(self):
        return self

    def __neg__(self):
        return self.times(-1)

    def __add__(self, other):
        assert self.n == other.n and self.m == other.m
        res = [[0.] * self.m for _ in range(self.n)]
        for i in range(self.n):
            res_i, self_i, other_i = res[i], self[i], other[i]
            for j in range(self.m):
                res_i[j] = (self_i[j] + other_i[j])
        return Matrix(self.n, self.m, res)

    def __sub__(self, other):
        assert self.n == other.n and self.m == other.m
        res = [[0.] * self.m for _ in range(self.n)]
        for i in range(self.n):
            res_i, self_i, other_i = res[i], self[i], other[i]
            for j in range(self.m):
                res_i[j] = (self_i[j] - other_i[j])
        return Matrix(self.n, self.m, res)

    def __mul__(self, other):
        if other.__class__ == Matrix:
            assert self.m == other.n
            res = [[0.] * other.m for _ in range(self.n)]
            for i in range(self.n):
                res_i, self_i = res[i], self[i]
                for k in range(self.m):
                    self_ik, other_k = self_i[k], other[k]
                    for j in range(other.m):
                        res_i[j] += self_ik * other_k[j]
            return Matrix(self.n, other.m, res)
        elif other.__class__ == Vector:
            assert self.m == other.n
            res = [0.] * self.n
            for i in range(self.n):
                for j in range(self.m):
                    res[i] += self[i][j] * other[j]
            return Vector(self.n, res)
        else:
            return self.times(other)

    def __rmul__(self, other):
        if other.__class__ == Vector: #転置をとってかける
            assert other.n == self.n
            res = [0.] * self.m
            for i in range(self.m):
                for j in range(self.n):
                    res[i] += self[j][i] * other[j]
            return Vector(self.m, res)
        return self.times(other)

    def __pow__(self, k):
        assert self.is_square()
        tmp = Matrix(self.n, self.n, self.mat)
        res = Matrix.id(self.n)
        while k:
            if k & 1:
                res *= tmp
            tmp *= tmp
            k >>= 1
        return res

    def determinant(self):
        assert self.is_square()
        res = 1
        tmp = Matrix(self.n, self.n, self.mat)
        for j in range(self.n):
            if is_zero(tmp[j][j]):
                for i in range(j + 1, self.n):
                    if not is_zero(tmp[i][j]):
                        break
                else:
                    return 0
                tmp.mat[j], tmp.mat[i] = tmp.mat[i], tmp.mat[j]
                res *= -1
            tmp_j = tmp[j]
            inv = 1 / tmp_j[j]
            for i in range(j + 1, self.n):
                tmp_i = tmp[i]
                c = -inv * tmp_i[j]
                for k in range(self.n):
                    tmp_i[k] += c * tmp_j[k]
        for i in range(self.n):
            res *= tmp[i][i]
        return res

    def inverse(self):
        assert self.is_square()
        res = Matrix.id(self.n)
        tmp = Matrix(self.n, self.n, self.mat)
        for j in range(self.n):
            if is_zero(tmp[j][j]):
                for i in range(j + 1, self.n):
                    if not is_zero(tmp[i][j]):
                        break
                else:
                    return -1
                tmp.mat[j], tmp.mat[i] = tmp.mat[i], tmp.mat[j]
                res.mat[j], res.mat[i] = res.mat[i], res.mat[j]
            tmp_j, res_j = tmp[j], res[j]
            inv = 1 / tmp_j[j]
            for k in range(self.n):
                tmp_j[k] *= inv
                res_j[k] *= inv
            for i in range(self.n):
                if i == j: continue
                c = tmp[i][j]
                tmp_i, res_i = tmp[i], res[i]
                for k in range(self.n):
                    tmp_i[k] -= tmp_j[k] * c
                    res_i[k] -= res_j[k] * c
        return res

    def linear_equations(self, vec):
        assert self.n == len(vec)
        aug = [self[i] + [vec[i]] for i in range(self.n)]
        rank = 0
        p = []
        q = []
        for j in range(self.m + 1):
            for i in range(rank, self.n):
                if not is_zero(aug[i][j]): break
            else:
                q.append(j)
                continue
            if j == self.m: return -1, [], []
            p.append(j)
            aug[rank], aug[i] = aug[i], aug[rank]
            inv = 1 / aug[rank][j]
            aug_rank = aug[rank]
            for k in range(self.m + 1):
                aug_rank[k] *= inv
            for i in range(self.n):
                if i == rank: continue
                aug_i = aug[i]
                c = -aug_i[j]
                for k in range(self.m + 1):
                    aug_i[k] += c * aug_rank[k]
            rank += 1
        dim = self.m - rank
        sol = [0.] * self.m
        for i in range(rank):
            sol[p[i]] = aug[i][-1]
        vecs = [[0.] * self.m for _ in range(dim)]
        for i in range(dim):
            vecs[i][q[i]] = 1.
        for i in range(dim):
            vecs_i = vecs[i]
            for j in range(rank):
                vecs_i[p[j]] = -aug[j][q[i]]
        return dim, sol, vecs

    def transpose(self):
        res = [[0.] * self.n for _ in range(self.m)]
        for i in range(self.n):
            for j in range(self.m):
                res[j][i] = self[i][j]
        return Matrix(self.m, self.n, res)

    def column(self, k):
        res = [0.] * self.n
        for i in range(self.n):
            res[i] = self[i][k]
        return Vector(self.n, res)

    def row(self, k):
        return Vector(self.m, self[k])

    def is_regular(self):
        if is_zero(self.determinant()):
            return False
        return True

class Vector():
    def __init__(self, n, vec=None):
        self.n = n
        self.vec = [0.] * n
        if vec:
            assert len(vec) == n
            self.vec = vec.copy()

    def __getitem__(self, key):
        if isinstance(key, slice):
            return self.vec[key]
        else:
            assert key >= 0
            return self.vec[key]

    def __setitem__(self, key, val):
        assert key >= 0
        self.vec[key] = val

    def __len__(self):
        return len(self.vec)

    def __str__(self):
        return ' '.join(map(str, self.vec))
    
    def __pos__(self):
        return self

    def __neg__(self):
        return self.times(-1)

    def __add__(self, other):
        assert self.n == other.n
        res = [0.] * self.n
        for i in range(self.n):
            res[i] = self[i] + other[i]
        return Vector(self.n, res)

    def __sub__(self, other):
        assert self.n == other.n
        res = [0.] * self.n
        for i in range(self.n):
            res[i] = self[i] - other[i]
        return Vector(self.n, res)

    def __mul__(self, other): 
        if other.__class__ == Vector: #内積
            assert self.n == other.n
            res = 0
            for i in range(self.n):
                res += self[i] * other[i]
            return res
        else:
            return self.times(other)

    def __rmul__(self, other):
        return self.times(other)

    def times(self, k):
        res = [0.] * self.n
        for i in range(self.n):
            res[i] = self[i] * k
        return Vector(self.n, res)

class SimplexMethod():
    # Ax <= b, x >= 0 の条件下で c * x を最大化する
    def __init__(self, c, A, b):
        self.n = len(c) + len(b) #変数の個数（スラック変数を含む）もとの変数は n - m 個
        self.m = len(b) #制約の個数
        self.A = A
        self.b = b
        self.c_B = Vector(self.m)
        self.c_N = Vector(self.n - self.m, c.vec)
        self.B = Matrix.id(self.m)
        self.N = Matrix(self.m, self.n - self.m, A.mat)
        self.B_idx = list(range(self.n - self.m, self.n))
        self.N_idx = list(range(self.n - self.m))
        self.xB = None
        self.x_opt = None

    def simplex(self):
        while True:
            invB = self.B.inverse()
            b_ = invB * self.b
            c_N_ = self.c_N - self.N.transpose() * (invB.transpose() * self.c_B)
            for i in range(self.n - self.m):
                if c_N_[i] > 0:
                    k = i # c_N_ が正になる一番最初の添字
                    break
            else:
                self.xB = b_
                return 1 # c_N_ <= 0 のとき最適解となる
            a_k = invB * self.N.column(k)
            for i in range(self.m):
                if a_k[i] > 0:
                    break
            else:
                return 0 # 非有界
            theta_min = float('inf')
            for i in range(self.m):
                if a_k[i] > 0:
                    theta = b_[i] / a_k[i]
                    if theta_min > theta:
                        theta_min = theta 
                        l = i # theta が最小になる一番最初の添字
            self.swapNB(l, k)

    def swapNB(self, l, k): # 基底変数 l と非基底変数 k を入れ替える
        self.N_idx[k], self.B_idx[l] = self.B_idx[l], self.N_idx[k]
        self.c_N[k], self.c_B[l] = self.c_B[l], self.c_N[k]
        N_k = self.N.column(k)
        B_l = self.B.column(l)
        for i in range(self.m):
            self.N[i][k] = B_l[i]
            self.B[i][l] = N_k[i]

    def find_feasible_solution(self):
        b_min = min(self.b)
        if b_min >= 0: return 1
        c_aux = Vector(self.n - self.m + 1, [-1] + [0] * (self.n - self.m))
        A_aux_mat = [[0.] * (self.n - self.m + 1) for _ in range(self.m)]
        for i in range(self.m):
            A_aux_mat[i][0] = -1
            for j in range(self.n - self.m):
                A_aux_mat[i][j + 1] = self.A[i][j]
        A_aux = Matrix(self.m, self.n - self.m + 1, A_aux_mat)
        S_aux = SimplexMethod(c_aux, A_aux, b)
        k = self.b.vec.index(b_min)
        S_aux.swapNB(k, 0)
        assert S_aux.simplex() # 非有界なことはあり得る？
        z_aux = S_aux.c_B * S_aux.xB
        if z_aux < 0:
            return 0 # 実行可能解がない
        if not 0 in S_aux.N_idx: # 退化して 0 が基底変数となっていた場合の処理は？
            r = S_aux.B_idx.index(0)
            for i in range(self.n - self.m + 1):
                S_aux.swapNB(r, i)
                if S_aux.B.is_regular():
                    break
                S_aux.swapNB(r, i)
        N_idx_aux = [S_aux.N_idx[i] - 1 for i in range(self.n - self.m + 1) if S_aux.N_idx[i] != 0]
        B_idx_aux = [S_aux.B_idx[i] - 1 for i in range(self.m)]
        idx_aux = N_idx_aux + B_idx_aux
        c_NB = self.c_N.vec + self.c_B.vec
        c_NB_sorted = [c_NB[idx_aux[i]] for i in range(self.n)]
        self.c_N = Vector(self.n - self.m, c_NB_sorted[:self.n - self.m])
        self.c_B = Vector(self.m, c_NB_sorted[self.n - self.m:])
        NB_mat = [[0.] * self.n for _ in range(self.m)]
        for i in range(self.m):
            for j in range(self.n - self.m):
                NB_mat[i][j] = self.N[i][j]
            for j in range(self.m):
                NB_mat[i][j + self.n - self.m] = self.B[i][j]
        NB_mat_sorted = [[0.] * self.n for _ in range(self.m)]
        for i in range(self.m):
            for j in range(self.n):
                NB_mat_sorted[i][j] = NB_mat[i][idx_aux[j]]
        N_mat = [[0.] * (self.n - self.m) for _ in range(self.m)]
        B_mat = [[0.] * self.m for _ in range(self.m)]
        for i in range(self.m):
            for j in range(self.n - self.m):
                N_mat[i][j] = NB_mat_sorted[i][j]
            for j in range(self.m):
                B_mat[i][j] = NB_mat_sorted[i][j + self.n - self.m]
        self.N = Matrix(self.m, self.n - self.m, N_mat)
        self.B = Matrix(self.m, self.m, B_mat)
        self.N_idx = N_idx_aux
        self.B_idx = B_idx_aux
        return 1

    def solve(self):
        if not self.find_feasible_solution(): raise Exception('Infeasible')
        if not self.simplex(): raise Exception('Unbounded')
        z = self.c_B * self.xB
        x = Vector(self.n, self.xB.vec + [0] * (self.n - self.m))
        x_sorted = [0.] * (self.n)
        idx = self.B_idx + self.N_idx
        for i in range(self.n):
            x_sorted[idx[i]] = x[i]
        self.x_opt = Vector(self.n - self.m, x_sorted[:self.n - self.m])
        return z, self.x_opt

C, D = map(int, input().split())

c = Vector(2, [1000, 2000])
A = Matrix(2, 2, [[3 / 4, 2 / 7], [1 / 4, 5 / 7]])
b = Vector(2, [C, D])

S = SimplexMethod(c, A, b)

z, x = S.solve()

print(z)
yukicoder

結果

テストケース

ソースコード
