#!/usr/bin/env python3

# author vas 17226

import numpy as np
from fractions import Fraction
import matplotlib.pyplot as plt
from copy import *
import math

def sgntostr(sgn):
    if sgn < 0:
        return "<"
    elif sgn > 0:
        return ">"
    else:
        return "="

class LPT:
    def __init__(self, kind, c, A, b, lim_signs, free_nrs, letter = 'x'):
        if kind != 'min' and kind != 'max':
            raise SystemError("Invalid kind")
        self.kind = kind
        self.n = A.shape[1]
        self.m = A.shape[0]
        self.c = c
        self.A = A
        self.b = b
        self.lim_signs = lim_signs
        self.free_nrs = list(map(lambda idx : idx - 1, free_nrs))
        self.letter = letter
        self.__xs_strs = [str(i + 1) for i in range(0, self.n)]

    def zero_var(self, idx):
        self.n -= 1
        self.c = np.delete(self.c, idx)
        self.A = np.delete(self.A, idx, axis=1)
        if idx in self.free_nrs:
            self.free_nrs = self.free_nrs.remove(idx)
        del self.__xs_strs[idx]

    def xs_strs(self, i):
        return self.letter + self.__xs_strs[i]

    def xs_comb_str(self, coefs):
        res =''
        for i in range(0, self.n):
            sign = ' +'
            if coefs[i] < 0:
                sign = ' -'
            res += sign + str(abs(coefs[i])) + self.xs_strs(i)
        return res

    def __str__(self):
        res = ''
        res += 'N=' + str(self.n) + " ; M=" + str(self.m) + '\n'
        res += self.xs_comb_str(self.c) + ' -> ' + self.kind + '\n'
        for rnum in range(0, self.A.shape[0]):
            res += self.xs_comb_str(self.A[rnum])
            sign = ''
            if self.lim_signs[rnum] == -1:
                sign = ' <= '
            elif self.lim_signs[rnum] == 0:
                sign = ' == '
            elif self.lim_signs[rnum] == 1:
                sign = ' >= '
            res += sign
            res += str(self.b[rnum])
            res += '\n'

        for i in range(0, self.n):
            isfree = ''
            if i in self.free_nrs:
                isfree = self.xs_strs(i) + ' free'
            else:
                isfree = self.xs_strs(i) + ' >= 0'
            res += isfree + '\n'

        return res

    def dual_classic_straight(self):
        res = copy(self)
        for rnum in range(0, res.m):
            if (res.kind == 'min' and res.lim_signs[rnum] == -1) or (res.kind == 'max' and res.lim_signs[rnum] == 1):
                res.A[rnum] *= -1
                res.b[rnum] *= -1
                res.lim_signs[rnum] *= -1

        return res

    def dual(self):
        res = self.dual_classic_straight()

        old_free_nrs = res.free_nrs.copy()

        res.free_nrs = []
        for i in range(0, res.m):
            if res.lim_signs[i] == 0:
                res.free_nrs += [i]

        res.lim_signs = [0 for i in range(0, res.n)]
        for i in range(0, res.n):
            if not (i in old_free_nrs):
                if res.kind == 'min':
                    res.lim_signs[i] = -1
                elif res.kind == 'max':
                    res.lim_signs[i] = 1
                else:
                    raise SystemError('haha lol')

        res.A = res.A.transpose()
        res.n, res.m = res.m, res.n
        res.c, res.b = res.b, res.c
        if res.kind == 'min':
            res.kind = 'max'
        else:
            res.kind = 'min'

        if res.letter == 'x':
            res.letter = 'y'
        else:
            res.letter = 'x'

        res.__xs_strs = [str(i + 1) for i in range(0, res.n)]

        return res

    def get_lim_signs(self, vec):
        lims = self.A.dot(vec)
        res = [0 for i in range(0, self.m)]
        for i in range(0, self.m):
            if lims[i] > self.b[i]:
                res[i] = 1
            elif lims[i] < self.b[i]:
                res[i] = -1
            else:
                res[i] = 0
        return lims, res


    def is_feasible(self, vec):
        for i in range(0, self.n):
            if (not (i in self.free_nrs)) and vec[i] < 0:
                return False, self.xs_strs(i) + ' < 0 (J1 failed)', 0

        lims, lim_signs_vec = self.get_lim_signs(vec)
        for i in range(0, self.m):
            if self.lim_signs[i] != lim_signs_vec[i] and (lim_signs_vec[i] != 0):
                err = 'For limitation #%s (counting from 1): %s [%s, should be %s] %s' \
                      % ((i+1), lims[i], lim_signs_vec[i], self.lim_signs[i], self.b[i])
                return False, err, 0

        return True, ','.join(map(sgntostr, lim_signs_vec)), self.c.dot(vec)

    def for_vec(self, vec):
        if len(vec) != self.n:
            raise SystemError('lenvec != self.n')
        msg = ''
        msg += 'For vector ' + str(vec) + ':\n'
        is_fsb, msg_fsb, obj_val = self.is_feasible(vec)
        msg += 'Feasibility: ' + str((is_fsb, msg_fsb, obj_val))
        msg += '\n'
        print(msg)
        return is_fsb, obj_val


def find_feasible_basis(A, b):
    At = A.transpose()
    ncols = At.shape[0]
    any_found = False
    for i in range(0, ncols):
        for j in range(i + 1, ncols):
            for k in range(j + 1, ncols):
                basis = np.array([At[i], At[j], At[k]])
                try:
                    det = np.linalg.det(basis)
                    print('det: %s' % det)
                    print('i j k: %s %s %s' % (i + 1, j + 1, k + 1))
                    print('B=\n' + str(basis.transpose()))
                    print('B^-1 * b=' + str(basis_inv.dot(b)))
                    basis_inv = np.linalg.inv(basis)
                    if abs(det) > 1e-2 and np.all(np.greater_equal(basis_inv.dot(b), np.zeros(b.shape))):
                        any_found = True
                        print('det = ' + str(det))
                    print('')
                except:
                    continue

    if not any_found:
        print('No feasible basis found')


# sols' rows are solutions
def is_solutions_optimal():
    xs_to_check = np.array([
        [0, 5./2, 1./2, -1./2],
        [1., 0, 0, -2.],
        [0, 0., 0, 1]
    ])

    kind = 'min'
    c = np.array([1., 5, -3, -5])
    A = np.array([
        [-2., -1, 4, 1],
        [1, -1, 1, 2],
        [3, 2, -1, -1]
    ])
    b = np.array([1., -3, 5])
    lim_signs = np.array([-1, -1, 0])
    free_nrs = [1, 4]

    ############

    straight = LPT(kind, c, A, b, lim_signs, free_nrs)
    dual = straight.dual()
    print('Straight: ', straight, '\n', 'Dual: ', dual, '')

    min_fsb_x = None
    min_func_obj = 10000000
    for x in xs_to_check:
        is_fsb, func_obj = straight.for_vec(x)
        if is_fsb and func_obj < min_func_obj:
            min_func_obj = func_obj
            min_fsb_x = x

    if min_fsb_x is None:
        print('No feasible X found!')
        return

    print('\nSolving for x=%s,  f(x)=%s' % (min_fsb_x, min_func_obj))

    # True, ','.join(map(sgntostr, lim_signs_vec)), self.c.dot(vec)
    lim_signs = straight.get_lim_signs(min_fsb_x)[1]
    print('lim_signs (Написать, что прямые ограничения выполняются со знаком): %s' % lim_signs)

    y_test = [0. for i in range(0, dual.n)]

    for i in range(0, straight.n):
        if abs(min_fsb_x[i]) > 1e-2: # then c_j = y*A_j
            print('По 2ой теореме выполняется равенство в двойственном ограничении #%s' % (i+1))
            dual.lim_signs[i] = 0

    for i in range(0, straight.m):
        if lim_signs[i] != 0:
            print('По 2ой теореме, так как прямое ограничение выполняется как неравенство, то y%s = 0' % (i+1))
            dual.zero_var(i)

    print('\nNew dual: ', dual)

    n_equalities = 0
    equality_A = np.zeros((dual.n, dual.n))
    equality_b = np.array((dual.n, 1))
    for i in range(0, dual.m):
        if dual.lim_signs[i] == 0:
            equality_A[n_equalities] = dual.A[i]
            equality_b[n_equalities] = dual.b[i]
            n_equalities += 1

    if n_equalities == dual.n:
        print('Можно решить систему')
        ys = np.linalg.solve(equality_A, equality_b)
        print('Игреки: ' + str(ys))


def dual_geom_plot():
    kind = 'min'
    c = np.array([1., 1.])
    A = np.array([
        [3., 2.],
        [2., 2.],
        [1., -1.],
        [4., -1.]
    ])
    b = np.array([6., 3., -1., -2.])
    lim_signs = np.array([0, 1, 1, 1])
    free_nrs = [2]

    ######

    straight = LPT(kind, c, A, b, lim_signs, free_nrs)
    dual = straight.dual()
    print('Straight: ', straight, '\n', 'Dual: ', dual, '')

    if dual.n != 2:
        raise ValueError("dual.n != 2 !!!")

    t = np.arange(-5., 5., 0.1)
    # y1 - horizontal , y2 vertival
    patches = []
    plt.axis('equal')
    for i in range(0, dual.m):
        a = dual.A[i, 0]
        b = dual.A[i, 1]
        c = dual.b[i]
        import matplotlib.patches as mpatches
        color=np.random.rand(3,)
        patch = mpatches.Patch(label='(%s) y1 + (%s) y2 ?? %s' % (a, b, c), color=color)
        patches += [patch]
        plt.plot(t, (c - a*t)/b, color=color)
    plt.legend(handles=patches)
    plt.grid()
    plt.show()


def example_1():
    kind = 'max'
    c = np.array([6, 3, -1, -2])
    A = np.array([[3, 2, 1, 4],
                  [2, 2, -1, -1]])
    b = np.array([0, 1])
    lim_signs = np.array([-1, 0])
    free_nrs = [1]

    lpt = LPT(kind, c, A, b, lim_signs, free_nrs)

    print('STRAIGHT: ' + str(lpt))

    lpt.for_vec([0, 2.5, 0.5, -0.5])
    lpt.for_vec([1, 0, 0, -2])
    lpt.for_vec([0, 0, 0, 1])

    lpt2 = lpt.dual()
    print('DUAL: ' + str(lpt2))



class STable:
    def __init__(self, T, initial_sigma, initial_omega):
        self.T = T.copy().astype('object')
        for i in range(0,self.T.shape[0]):
            for j in range(0, self.T.shape[1]):
                self.T[i,j]=Fraction(int(self.T[i,j]))

        self.sigma = initial_sigma
        self.omega = initial_omega
        if not self.straight_dop():
            raise SystemError('not straight dop')

    def straight_dop(self):
        col = np.array(self.T[1:,0])
        return np.all(np.greater_equal(col, np.zeros(col.shape)))

    def dual_dop(self):
        row = np.array(self.T[0,1:])
        return np.all(np.greater_equal(row, np.zeros(row.shape)))

    def doswitch(self, rnum, snum):
        for i in range(0, self.T.shape[0]):
            if i != rnum:
                coef = self.T[i,snum] / self.T[rnum,snum]
                self.T[i] -= self.T[rnum] * coef

        self.T[rnum] /= self.T[rnum,snum]
        self.sigma[rnum] = self.omega[snum]

    def dancig_choose(self):
        snum = 1
        min_val = 100000000.
        min_val_snum = 1
        while snum < self.T.shape[1]:
            if self.T[0,snum] < min_val and self.T[0, snum] < 0:
                min_val = self.omega[snum]
                min_val_snum = snum
            snum += 1

        rnum = 1
        min_coef = 100000000.
        min_coef_rnum = 1
        while rnum < self.T.shape[0]:
            coef = self.T[rnum,0] / self.T[rnum, min_val_snum]
            if coef < min_coef:
                min_coef = coef
                min_coef_rnum = rnum
            rnum += 1

        return min_coef_rnum, min_val_snum

    def blend_choose(self):
        snum = 1
        min_omega = 1000000
        min_omega_snum = 1
        while snum < self.T.shape[1]:
            if self.omega[snum] < min_omega and self.T[0,snum] < 0:
                min_omega = self.omega[snum]
                min_omega_snum = snum
            snum += 1

        rnum = 1
        min_sigma = 1000000
        min_sigma_rnum = 1
        while rnum < self.T.shape[0]:
            if self.sigma[rnum] < min_sigma and self.T[rnum,min_omega_snum] > 0:
                min_sigma = self.sigma[rnum]
                min_sigma_rnum = rnum
            rnum += 1

        return min_sigma_rnum, min_omega_snum

    def bdr(self):
        res = [0.0 for i in range(1,len(self.omega))]
        for rnum in range(1, self.T.shape[0]):
            res[self.sigma[rnum]] = self.T[rnum, 0]
        return res

    def smethod(self):
        while True:
            print('-----------------------------------------------------------')
            print('omega: %s\nsigma: %s\n' % (self.omega[1:], self.sigma[1:]))
            table_str = ''
            for r in self.T:
                for v in r:
                    if v.denominator == 1:
                        table_str += '%6s' % str(v.numerator)
                    else:
                        table_str += '%6s' % (str(v.numerator) + '/' + str(v.denominator))
                    table_str += ' '
                table_str += '\n'
            print('table: \n%s' % table_str)

            if self.dual_dop():

                print('\nFound solution=%s at vector %s' % (-self.T[0,0], self.bdr()))
                return

            rnum, snum = self.blend_choose()
            # rnum, snum = self.dancig_choose()
            print('\nUsing [row=%s, col=%s] (numerating from 0)' % (rnum, snum))

            if np.all(np.less_equal(self.T[1:,snum], np.zeros(self.T[1:,snum].shape))):
                print('\nTheres no solution!!')
                return

            self.doswitch(rnum, snum)



def example_2():
    T = np.array([
        [-6, 0, 0, 0, -3, 1],
        [3, 0, 0, 1, 4, -2],
        [8, 1, 0, 0, 8, 1],
        [16, 0, 1, 0, 19, -2]
    ])

    # first value is omitted
    sigma = np.array([0,3,1,2])

    # first value is omitted
    omega = np.array([0,1,2,3,4,5])

    st = STable(T, sigma, omega)
    st.smethod()

def slayter(test_pt, phis):
    for phi in phis:
        if phi(test_pt) >= 0:
            return False

        return True

def find_active(test_pt, phis):
    active_idxs = []
    for i in range(0, len(phis)):
        phi = phis[i]
        if abs(phi(test_pt)) < 1e-4:
            active_idxs += [i]
    return active_idxs

def select_by_indices(arr, indices):
    res = []
    for i in indices:
        res += [arr[i]]
    return res

def example_3():

    check_points = [
        [2.,math.log(2.)],
        [math.e, 0.],
        [1.,0.]
    ]

    phis = [
        (lambda x: -math.log(x[0]) + x[1]),
        (lambda x: -x[0]),
        (lambda x: -x[1])
    ]

    slayter_pt = np.array([math.exp(2), 0.5])

    # производная f
    f_d = lambda x: [math.exp(x[0]+x[1]) + 2.*x[0],
                     math.exp(x[0]+x[1]) - 2.]

    phi_ds = [
        lambda x: [ #Phi1
            -1./x[0],
            1.
        ],
        lambda x: [ #Phi2
            -1.,
            0.
        ],
        lambda x: [
            0.,
            -1.
        ]
    ]

    print('Slayter: ', slayter(slayter_pt, phis))
    for xs in check_points:
        print('\ntesting point %s' % xs)

        active_idxs = find_active(xs, phis)
        print('active idxs: %s' % list(map(lambda idx : idx+1, active_idxs)))
        b = -1.*np.array(f_d(xs))
        A = np.column_stack(tuple(
            list(phi_i_d(xs) for phi_i_d in select_by_indices(phi_ds, active_idxs))
        ))
        print('A: \n%s\n b: %s' % (A, b))
        # lambdas = np.linalg.solve(A, b)
        # print('Lambdas: ' + str(lambdas))

def main():
    # A = np.array([
    #     [2., -1., 1., 1., 1.],
    #     [-4., 3., -1., -1., -3.],
    #     [3., 2., 3., 5., 0.]
    # ])
    # print(A.shape)
    # b= np.array([2., -4., 3.])
    # find_feasible_basis(A, b)


    # is_solutions_optimal()
    dual_geom_plot()

main()