6
\$\begingroup\$

In my Uni, my scientific professor asked me to make some researches about the extreme points of polyhedrals. And I did them. I found that there is still no code in public for searching extreme points for polyhedral with n dimensions (n - x's), but polyhedrons are everywhere (CV, game theories, etc.). I wrote a function for this task and made a python library (also there are matrix and array combination maker functions).

All I want is to make this code more optimal and compatible for all python versions (I have some troubles that some times happen when I install it by "pip install lin", but other times no). I want to make the life of people easier and make it more comfortable.

I am asking for you to test this function on your computer and write if you have any bugs, fails or thoughts on how to make it better. I am open to constructive criticism and non-constructive too (it will help to understand if somebody needs it or that is just a waste of time).

All the examples, instructions and code on my GitHub: https://github.com/r4ndompuff/polyhedral_set

import numpy as np
import itertools as it
import math
import re

def permutation(m,n):
    return math.factorial(n)/(math.factorial(n-m)*math.factorial(m))

def matrix_combinations(matr,n):
    timed = list(map(list, it.combinations(matr, n)))
    for i in range(n):
        timed[i][i][i] = np.asscalar(timed[i][i][i])
    all = np.array(list(timed))
    return all

def array_combinations(arr,n):
    timed = list(map(list, it.combinations(arr, n)))
    for i in range(n):
        timed[i][i] = np.asscalar(timed[i][i])
    all = np.array(list(timed))
    return all

def check_extreme(matr, arr, x, sym_comb, m):
    sym_comb = sym_comb.replace(']', '')
    sym_comb = sym_comb.replace('[', '')
    sym_comb = re.split("[ ,]", sym_comb)
    for i in range(m):
        td_answer = sum(matr[i]*x)
        if sym_comb[i] == '>':
            if td_answer <= arr[i]:
                return 0
        elif sym_comb[i] == '>=':
            if td_answer < arr[i]:
                return 0
        elif sym_comb[i] == '<':
            if td_answer >= arr[i]:
                return 0
        elif sym_comb[i] == '<=':
            if td_answer > arr[i]:
                return 0
        elif sym_comb[i] == '=':
            if td_answer != arr[i]:
                return 0
        elif sym_comb[i] == '!=':
            if td_answer == arr[i]:
                return 0
        else:
            return 0
    return 1

def extreme_points(m,n,A,b,sym_comb):
    # Input
    A = np.array(A).reshape(m,n)
    b = np.array(b).reshape(m,1)
    # Proccess
    ans_comb = np.zeros((1,n))
    arr_comb = array_combinations(b,n)
    matr_comb = matrix_combinations(A,n)
    for i in range(int(permutation(n,m))):
        if np.linalg.det(matr_comb[i]) != 0:
            x = np.linalg.solve(matr_comb[i],arr_comb[i])
            ans_comb = np.vstack([ans_comb,x])
    ans_comb = np.delete(ans_comb, (0), axis=0)
    j = 0
    for i in range(len(ans_comb)):
        if check_extreme(A, b, ans_comb[j], sym_comb, m):
            ans_comb = ans_comb
            j = j + 1
        else:
            ans_comb = np.delete(ans_comb, (j), axis=0)
    # Output
    return ans_comb

And I am uploading some more tests. https://imgur.com/mjweDyy

\$\endgroup\$
  • \$\begingroup\$ Can you provide some more information on what exactly is going on, from a math standpoint? Is this actually a polyhedron, or a polytope? In how many dimensions? By "extreme", what do you mean? Euclidean norm from (the origin, an arbitrary point)? \$\endgroup\$ – Reinderien Apr 21 at 23:21
  • \$\begingroup\$ "there is still no code in public for searching extreme points for polyhedral with n dimensions" - I can nearly guarantee that that isn't the case. \$\endgroup\$ – Reinderien Apr 21 at 23:22
  • 2
    \$\begingroup\$ @Reinderien I am still in the process of deciphering the question. I have a math background, and I share the mother tongue with OP. My impression is that by extremal points OP means the vertices of a simplex where a certain linear form (defined by A, b, and the condition) achieves an extremum. I could be wrong. \$\endgroup\$ – vnp Apr 22 at 0:32
  • \$\begingroup\$ @vnp That's kind of what I guessed, and if that's the case, linear programming is a quite well-established field already - with some stuff built right into scipy. \$\endgroup\$ – Reinderien Apr 22 at 0:33
  • \$\begingroup\$ @Reinderien Agreed. Still deciphering. \$\endgroup\$ – vnp Apr 22 at 0:34
5
\$\begingroup\$

I created a rudimentary pull request to your GitHub repo. I won't show all of the content here except for the main file:

import numpy as np
import itertools as it
from math import factorial
import re


def permutation(m, n):
    return factorial(n) / (factorial(n - m) * factorial(m))


def matrix_combinations(matr, n):
    timed = list(map(list, it.combinations(matr, n)))
    return np.array(list(timed))


def array_combinations(arr, n):
    timed = list(map(list, it.combinations(arr, n)))
    return np.array(list(timed))


def check_extreme(matr, arr, x, sym_comb, m):
    sym_comb = sym_comb.replace(']', '')
    sym_comb = sym_comb.replace('[', '')
    sym_comb = re.split("[ ,]", sym_comb)
    for i in range(int(m)):
        td_answer = sum(matr[i] * x)
        if sym_comb[i] == '>':
            if td_answer <= arr[i]:
                return 0
        elif sym_comb[i] == '>=':
            if td_answer < arr[i]:
                return 0
        elif sym_comb[i] == '<':
            if td_answer >= arr[i]:
                return 0
        elif sym_comb[i] == '<=':
            if td_answer > arr[i]:
                return 0
        elif sym_comb[i] == '=':
            if td_answer != arr[i]:
                return 0
        elif sym_comb[i] == '!=':
            if td_answer == arr[i]:
                return 0
        else:
            return 0
    return 1


def extreme_points(A, b, sym_comb):
    # Input
    A = np.array(A)
    b = np.array(b)
    m, n = A.shape
    # Process
    ans_comb = np.zeros((1, n))
    arr_comb = array_combinations(b, n)
    matr_comb = matrix_combinations(A, n)
    for i in range(int(permutation(n, m))):
        if np.linalg.det(matr_comb[i]) != 0:
            x = np.linalg.solve(np.array(matr_comb[i], dtype='float'),
                                np.array(arr_comb[i], dtype='float'))
            ans_comb = np.vstack([ans_comb, x])
    ans_comb = np.delete(ans_comb, 0, axis=0)
    j = 0
    for i in range(len(ans_comb)):
        if check_extreme(A, b, ans_comb[j], sym_comb, m):
            ans_comb = ans_comb
            j += 1
        else:
            ans_comb = np.delete(ans_comb, j, axis=0)
    # Output
    return ans_comb

Notable changes:

  • Do a direct import of factorial
  • Don't call asscalar, since it's both unneeded and deprecated
  • Don't call a variable all, since that shadows a Python built-in
  • Don't need to explicitly pass array dimensions, nor do you need to reshape the arrays
  • Drop redundant parens around some expressions
  • Use += where applicable
  • Fix up almost all PEP8 issues, except for your capital letter A, which is fine in context

This doesn't solve the bigger issue that you should replace 99% of this with a call to scipy. I'll do that separately (I suspect that @vnp is, as well).

\$\endgroup\$
  • 1
    \$\begingroup\$ Thank you! In meantime I made it work with any NumPy version now (about loop error you report). \$\endgroup\$ – Andrey Lovyagin Apr 22 at 7:57
  • \$\begingroup\$ Sorry, but I can't find the same function in scipy. Can you provide me the name of function? \$\endgroup\$ – Andrey Lovyagin Apr 22 at 8:17
  • \$\begingroup\$ @AndreyLovyagin docs.scipy.org/doc/scipy/reference/generated/… \$\endgroup\$ – Reinderien Apr 22 at 16:18
  • \$\begingroup\$ Thank you very much for all that you have done! Also, I checked scipy.optimize.linprog() function and I can't really understand how it can replace much of the code. I will learn and read more about it, but on first sight, it looks like it can't solve non-equational system (I mean not '=', but any sign, e.g. '>='), that is why I can't see how I can find extreme points by it. \$\endgroup\$ – Andrey Lovyagin Apr 22 at 20:39
  • \$\begingroup\$ Read about the A_ub, b_ub arguments. \$\endgroup\$ – Reinderien Apr 22 at 20:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.