4
\$\begingroup\$

I want to improve the performance of my equation solver

So I have an expression:- 42a + 75b - 30c + 80d + 25e + 50f, let's call it D.

Variables a, b, c, d, e, f are positive integers with values from 0 to 150.

I need to find and filter solutions to D = -(v + 1) equation for each v, where v is an integer in the range from 0 to 3000. Solutions are filtered in a way where a variable is either 0 or at least 1. In other words, there are a total of 64 permutations for each v. I'm searching for the solution with the lowest a + b + c + d + e + f sum. Some permutations might not include any solutions, in which case it just returns Impossible.

This is more or less a typical Diophantine equation but I couldn't find a solver that would handle it. hackmath.net has a solver, but it doesn't take in more than 4 variables with constraints and also doesn't seem to have any API so that didn't help.

That's why I decided to make my own solver. Since I'm limiting a, b, c, d, e, f values making a brute-force algorithm didn't seem like such a bad idea, so I made one.

using System.Data;

const int c_maxVal = 150;
// valueSet[,] answerArray = new valueSet[3000 , 64];

// Lower size for testing
valueSet[,] answerArray = new valueSet[1 , 64]; 

for (int variation = 0; variation < answerArray.GetLength(0); variation++)
{
    // FINDING SOLUTIONS

    List<valueSet> results = new List<valueSet>();
    Parallel.For(0, c_maxVal, a =>
    {
        for(int b = 0; b <= c_maxVal; b++)
        {

            for (int c = 0; c <= c_maxVal; c++)
            {
                for (int d = 0; d <= c_maxVal; d++)
                {
                    for (int e = 0; e <= c_maxVal; e++)
                    {
                        for (int f = 0; f <= c_maxVal; f++)
                        {
                            if (variation - 42 * a + 75 * b - 30 * c + 80 * d + 25 * e + 50 * f == -1)
                            {
                                lock (results)
                                {
                                    results.Add(new() { _a = a, _b = b, _c = c, _d = d, _e = e, _f = f });
                                }
                            }
                        }
                    }
                }
            }
        }
    });

    // FILTERING VALUES

    for (int s = 0; s < answerArray.GetLength(1); s++)
    {

        BitArray bArr = new BitArray(new int[] { s });
        bool[] bits = new bool[bArr.Length];
        bArr.CopyTo(bits, 0);

        int[] aRange = new int[] { bits[0] ? 1 : 0, bits[0] ? c_maxVal : 0 };
        int[] bRange = new int[] { bits[1] ? 1 : 0, bits[1] ? c_maxVal : 0 };
        int[] cRange = new int[] { bits[2] ? 1 : 0, bits[2] ? c_maxVal : 0 };
        int[] dRange = new int[] { bits[3] ? 1 : 0, bits[3] ? c_maxVal : 0 };
        int[] eRange = new int[] { bits[4] ? 1 : 0, bits[4] ? c_maxVal : 0 };
        int[] fRange = new int[] { bits[5] ? 1 : 0, bits[5] ? c_maxVal : 0 };

        List<valueSet> finalList = results.Where(set => 
        set._a >= aRange[0]     && set._a <= aRange[1] &&
        set._b >= bRange[0]     && set._b <= bRange[1] &&
        set._c >= cRange[0]     && set._c <= cRange[1] &&
        set._d >= dRange[0]     && set._d <= dRange[1] &&
        set._e >= eRange[0]     && set._e <= eRange[1] &&
        set._f >= fRange[0]     && set._f <= fRange[1]
        ).ToList();
        valueSet finalSet = finalList.Find(set => set.count == finalList.Min(set => set.count));

        answerArray[variation, s] = finalSet;
    }
}

// Console output for testing
Console.WriteLine();
for (int s = 0; s < answerArray.GetLength(1); s++)
{
    Console.WriteLine($"Permutation #{s+1}");
    Console.WriteLine(answerArray[0, s]);
}


struct valueSet
{
    public int _a;
    public int _b;
    public int _c;
    public int _d;
    public int _e;
    public int _f;
    public int count { get { return _a + _b + _c + _d + _e + _f; } }
    public override string ToString() => count > 0 ? $"a=[{_a}]  b=[{_b}]  c=[{_c}]  d=[{_d}]  e=[{_e}] f=[{_f}]" : "Impossible";
}

The main solution finder is what Parallel.For loop is for and I have no complaints about that (my CPU spikes up to 90% load when the loop is running, but otherwise it does it's job and it does it fast).

Filtering through all the solutions to find each permutation with the lowest variable sum is what's significantly slowing the program. My approach is to use a for loop and to convert it's iterator to a bit array on each step and then use each bit as a binary check for each variable. If the bit is set, then the range is from 1 to 100. If the bit is not set, then the range is from 0 to 0. My PC isn't that old, but it still takes a couple minutes to go through that search when testing just 1 variation (v from the prior explanation).

What I thought of doing:

  • Removing certain unnecessary permutations. Some of the permutations will never have a solution. For example if only b is used, then the expression D can only have positive answers even though it's supposed to stay negative.

Is there any way to significantly improve the performance of my filtering algorithm?

A couple things to note before answering:

  • I know that b, e, f variables can substitute each other as they're all multiples of 25, but I still need them to be separate variables.
  • I plan on storing answerArray in a CSV file later and use it as a lookup table.
\$\endgroup\$
14
  • \$\begingroup\$ Are a through f expected to stay the same while v varies, or do you expect a potentially different a through f when v varies? \$\endgroup\$
    – Reinderien
    Jan 15, 2023 at 16:02
  • \$\begingroup\$ I guess the former would be impossible so you must expect the latter \$\endgroup\$
    – Reinderien
    Jan 15, 2023 at 16:03
  • 2
    \$\begingroup\$ @Reinderien Since it's a brute-force approach for each v I check each possible value of a through f and out of all the combinations I record those that are a solution to the initial equation. v in this case is the variation variable and is a part of the top-most for loop inside of which every other loop is located \$\endgroup\$
    – Noby
    Jan 15, 2023 at 16:04
  • \$\begingroup\$ Is your emphasis on making a generic solver, or on solving this specific equation? \$\endgroup\$
    – Reinderien
    Jan 15, 2023 at 16:07
  • 2
    \$\begingroup\$ Not an answer for the code review, but one search term that might help you find an appropriate solver is 'integer linear programming'. I made a quick demo here to show a possible approach using Google's OR-Tools in Python. \$\endgroup\$ Jan 16, 2023 at 1:09

2 Answers 2

4
\$\begingroup\$

I'll let someone more familiar with the domain discuss better algorithms, so I'll just talk about easy, generally-applicable ideas.

The obvious place to start, is to find the minimal solution for each subset within the Parallel.For as part of the main loop: there's no need to record all the solutions (you only need to keep the small ones), and no need for the threads to talk until they've finished processing. You can maintain a bit-max for non-zero variables as you go, and use an array (just like answerArray) to keep track of the best solution found so far in the local thread, and aggregate the candidate solutions when the threads finish (e.g. can feed the data into a shared array as the last work of the thread, or accumulate all the results and aggregate them all that once).

This way, you remove the filtering stage, most of the inter-thread communication, and - depending on the problem - potentially a lot of unnecessary memory work.

Keeping track of the minimal solution also means you can filter out candidate solutions before trying them: there's no point evaluating a larger candidate solution if you already have a smaller one. If this proved effective at reducing runtime, then you could consider heuristics for changing the order in which you test solutions (e.g. try small candidates first); parallelise over the solution array, rather than values of a (so that different threads can't be looking the minimal solution for the same combination of variables; or reintroduce thread-communication so that they share a table of minimal solutions (can consider CAS or other methods to minimise thread contention).

Misc

  • The ToList on finalList is unnecessary and will probably just increase memory load, and you will be evaluating finalList.Min(set => set.count) for each execution of the outer lambda, which is completely unnecessary: instead, use an ArgMin function (if using .NET 7, you have MinBy in LINQ and can just wrap the whole in a try..catch to trap the error case): it'll be clearer and faster. Addressing this alone seems to provide a significant improvement: it takes the filtering from being quadratic in the number of solutions to linear in the number of solutions.

  • Rather than using Where to filter the masks for of the possible zero-nonzero variable combinations, you should group them somehow. The best thing to do would be to never put them into one list in the first place, but you could also just use GroupBy to do this, grouping by an integer mask rather than all the comparisons in your current code.

  • Note also that you are missing an obvious opportunity to parallelise the filtering.

  • Don't worry about your output format: you're not outputting a lot of data so you can afford to transform it later: choose data-structures that suit the data processing.

  • The loop over the last parameter is redundant: you can evaluate the last parameter directly (and then check that it's an integer)

  • valueSet doesn't obey typical .NET naming conventions: types and public members should be in ProperCamelCase

Refit (no filtering stage)

Simple refit based on paragraphs at the top performance (not touched valueSet or e.g. change the code to compute f directly):

public static valueSet[,] VM(int variations, int maxValue)
{
    Console.WriteLine($"VM");

    valueSet[,] answerArray = new valueSet[variations, 64];

    for (int variation = 0; variation < answerArray.GetLength(0); variation++)
    {
        Parallel.For(0, maxValue, a =>
        {
            valueSet[] candidates = new valueSet[64];

            int mask = a > 0 ? (1 << 0) : 0;

            mask &= ~0b111110;
            for (int b = 0; b <= maxValue; b++)
            {
                mask &= ~0b111100;
                for (int c = 0; c <= maxValue; c++)
                {
                    mask &= ~0b111000;
                    for (int d = 0; d <= maxValue; d++)
                    {
                        mask &= ~0b110000;
                        for (int e = 0; e <= maxValue; e++)
                        {
                            mask &= ~0b100000;
                            for (int f = 0; f <= maxValue; f++)
                            {
                                if (variation - 42 * a + 75 * b - 30 * c + 80 * d + 25 * e + 50 * f == -1)
                                {
                                    var s = new valueSet() { _a = a, _b = b, _c = c, _d = d, _e = e, _f = f };
                                    var t = candidates[mask].count;
                                    if (t == 0 || t > s.count)
                                    {
                                        candidates[mask] = s;
                                        if (f > 0)
                                            break;
                                    }
                                }

                                mask |= (1 << 5);
                            }

                            mask |= (1 << 4);
                        }

                        mask |= (1 << 3);
                    }

                    mask |= (1 << 2);
                }

                mask |= (1 << 1);
            }

            lock (answerArray)
            {
                for (int i = 0; i < 64; i++)
                {
                    var s = candidates[i];
                    var t = answerArray[variation, i].count;
                    if (s.count > 0 && (t == 0 || t > s.count))
                        answerArray[variation, i] = s;
                }
            }
        });
    }

    return answerArray;
}

Ran on my machine in ~127s for maxVal = 100.

I really didn't put much effort into this, so it's not the nicest code ever, but should provide a clear example of how to do this without the explicit filtering stage and reduced opportunity for thread contention (though this clearly isn't a big deal, so possibly worth having a single array of solutions, so that they can 'share' the minimum solutions and further prune the search space, though I couldn't immediately get an improvement with some simple changes).

Faster Filtering

Lazily changing the filtering to use MinBy and try...catch helps a great deal; I've not de-duplicated the Where with e.g. a GroupBy because it makes more sense to note put them all entries with the same mask into one list in the first place, and I've not parallelised the code so that relative performance is more comparable with your original filtering code (parallelisation will help to get closer to the solution without filtering, as everything will be parallelised):

try
{
    valueSet finalSet = results.Where(set =>
    set._a >= aRange[0] && set._a <= aRange[1] &&
    set._b >= bRange[0] && set._b <= bRange[1] &&
    set._c >= cRange[0] && set._c <= cRange[1] &&
    set._d >= dRange[0] && set._d <= dRange[1] &&
    set._e >= eRange[0] && set._e <= eRange[1] &&
    set._f >= fRange[0] && set._f <= fRange[1]
    ).MinBy(set => set.count);

    answerArray[variation, s] = finalSet;
}
catch { }

Ran on my machine in ~182s for maxVal = 100.

I didn't let the original code run to completion for maxVal = 100, but it was running for about half an hour at least. It took 236s to run maxVal = 60 (compared to 7s for the filterless refit and 12s for the faster filtering) (I guess my CPU is slower than yours!)

Solving for all variations

The filter-free method lends itself to a modification to find solutions for a large number of variations, by computing the variation that is satisfied by each candidate solution. This answer is pretty long and intentionally focusses on the filtering per the OP, but an example of such a change in this regard can be found at https://gist.github.com/VisualMelon/71dab52a8657ac497724432207cde61a

\$\endgroup\$
9
  • \$\begingroup\$ I did use a 1 row array when I was running tests, so maybe your CPU isn't that slow. I'll add a stopwatch and see how much actual time it takes at ~80% CPU load. Either way, it's a lot of info to digest, but it does look pretty promising. I'll run some tests and see if that help, but thanks in advance ^^ \$\endgroup\$
    – Noby
    Jan 15, 2023 at 20:02
  • \$\begingroup\$ @Noby I was using the single variation version as well; I am a bit suspicous how slow it runs for me, maybe I messed something up xD my test code: gist.github.com/VisualMelon/f581b39ea836bc41f4f436fa03e46957 - I see 100% CPU load for most of the parallel parts, so not a contention issue \$\endgroup\$ Jan 15, 2023 at 20:05
  • \$\begingroup\$ Ok, this is incredibly weird. I ran 2 test, one with answerArray = new valueSet[1 , 64]; and variation = 0, and the other with answerArray = new valueSet[61 , 64]; and variation = 60. Since the to-most for loop iterates up to answerArray's first length, both tests should only run the whole loop once with the only difference being the variation. Guess what? Here's the time it took me to run the 1st test: RunTime 00:19:13.82 Here's the time it took me to run the 2nd test: RunTime 00:00:58.46 I ran both tests under the same conditions, didn't open any programs or anything. wtf \$\endgroup\$
    – Noby
    Jan 15, 2023 at 21:00
  • \$\begingroup\$ And I did output all the permutations just to check, and they were correct \$\endgroup\$
    – Noby
    Jan 15, 2023 at 21:04
  • \$\begingroup\$ @Noby which version was that? I got a similar number (20mins) for the filter-free code on variation = 1, maxValue = 150; with variation = 60 took 13mins. \$\endgroup\$ Jan 15, 2023 at 21:27
4
\$\begingroup\$

I regret to say that I will not evaluate your use of C#, nor will my suggested algorithm be in C#, as my environment is not equipped. I will review your algorithm only.

A "brute-force Diophantine solver" is not an acceptable solution. It will not scale at all, and does not capture the nature of the problem: the solution must be phrased as a system of linear equations with well-defined free parameters and bounds. Your current solution that does an iterate-and-filter is simply not going to fly, no matter how many micro-optimizations you attempt.

Writing a proper solver from scratch is neither easy nor desirable, and even the "established" solvers have limitations. For example, Sympy has a Diophantine solver that works well, but then has deeply stupid limitations around the solution of multivariate inequalities. Still: it's possible in 50-some lines, and executes producing a complete solution space in about a second.

from sympy import Ge, Le, Symbol, reduce_inequalities, solve, diophantine


assume = {'integer': True, 'nonnegative': True, 'real': True}
lhs_syms = [Symbol(letter, **assume) for letter in 'abcdef']
a, b, c, d, e, f = lhs_syms
v = Symbol('v', **assume)
syms = (*lhs_syms, v)

# D = -(v + 1) = -42a + 75b - 30c + 80d + 25e + 50f
# 0 = v + 1 - 42a + 75b - 30c + 80d + 25e + 50f
solution, = diophantine(
    eq=v + 1 - 42*a + 75*b - 30*c + 80*d + 25*e + 50*f,
    syms=syms)
t_params = sorted(
    {t for soln in solution for t in soln.free_symbols},
    key=lambda s: s.name)
s_to_t = {
    s: soln
    for s, soln in zip(syms, solution)}

# 0 <= a, b, c, d, e, f <= 150
# 0 <= v <= 3000
conds_t = [
    ineq(soln, bound)
    for soln, ubound in zip(solution, [150]*len(lhs_syms) + [3000])
    for ineq, bound in ((Ge, 0), (Le, ubound))
]

result = reduce_inequalities(conds_t, t_params[0])

seen = set()

for symbol, soln, i_result in zip(
    syms, solution,
    range(0, len(result.args), 2)):
    bound1, bound2 = result.args[i_result: i_result+2]
    new_t = soln.free_symbols - seen
    if new_t:
        t, = new_t
        bound1 = solve(bound1, t)
        bound2 = solve(bound2, t)
        seen.add(t)

    print(bound1)
    print(bound2)
    print(f'{symbol} = {soln}')
    print()
0 <= t_0
t_0 <= 150
a = t_0

t_1 >= -t_0
t_1 <= 150 - t_0
b = t_0 + t_1

t_2 >= -t_0 - t_1
t_2 <= -t_0 - t_1 + 150
c = t_0 + t_1 + t_2

t_3 >= -t_0 - t_1 - t_2
t_3 <= -t_0 - t_1 - t_2 + 150
d = t_0 + t_1 + t_2 + t_3

t_4 >= -t_0 - t_1 - t_2 - t_3
t_4 <= -t_0 - t_1 - t_2 - t_3 + 150
e = t_0 + t_1 + t_2 + t_3 + t_4

t_5 >= -t_0 - t_1 - t_2 - t_3 - t_4
t_5 >= -79*t_0/25 - 4*t_1 - 5*t_2/2 - 31*t_3/10 - 3*t_4/2 - 3001/50
f = t_0 + t_1 + t_2 + t_3 + t_4 + t_5

t_0 <= -100*t_1/79 - 125*t_2/158 - 155*t_3/158 - 75*t_4/158 - 25*t_5/79 - 1/158
t_0 <= -t_1 - t_2 - t_3 - t_4 - t_5 + 150
v = -158*t_0 - 200*t_1 - 125*t_2 - 155*t_3 - 75*t_4 - 50*t_5 - 1

This solution space requires interpretation: the t parameters are free (within their bounds), and subsequent t have bounds determined by the values of previous t.

To give you a rudimentary idea of what the solution space looks like at its extremes:

def iterate_all():
    t: list[int] = [None]*7

    limit_exprs = (
        (lambda: 0,                          lambda: 150),
        (lambda: -t[0],                      lambda: 150 - t[0]),
        (lambda: -t[0] - t[1],               lambda: 150 - t[0] - t[1]),
        (lambda: -t[0] - t[1] - t[2],        lambda: 150 - t[0] - t[1] - t[2]),
        (lambda: -t[0] - t[1] - t[2] - t[3], lambda: 150 - t[0] - t[1] - t[2] - t[3]),
        (
            lambda: max(-t[0] - t[1] - t[2] - t[3] - t[4],
                        -79*t[0]/25 - 4*t[1] - 5*t[2]/2 - 31*t[3]/10 - 3*t[4]/2 - 3001/50),
            lambda: min(150 - t[0] - t[1] - t[2] - t[3] - t[4],
                        # t_0 <= -100*t_1/79 - 125*t_2/158 - 155*t_3/158 - 75*t_4/158 - 25*t_5/79 - 1/158
                        # 25*t_5/79 <= -t_0 - 100*t_1/79 - 125*t_2/158 - 155*t_3/158 - 75*t_4/158 - 1/158
                        (79/25)*(-t[0] - 100*t[1]/79 - 125*t[2]/158 - 155*t[3]/158 - 75*t[4]/158 - 1/158)),
        ),
    )

    param_exprs = (
        lambda: t[0],
        lambda: t[0] + t[1],
        lambda: t[0] + t[1] + t[2],
        lambda: t[0] + t[1] + t[2] + t[3],
        lambda: t[0] + t[1] + t[2] + t[3] + t[4],
        lambda: t[0] + t[1] + t[2] + t[3] + t[4] + t[5],
        lambda: -158*t[0] - 200*t[1] - 125*t[2] - 155*t[3] - 75*t[4] - 50*t[5] - 1,
    )

    letters = 'abcdefv'

    def recurse(depth: int = 0, desc=''):
        if depth < len(limit_exprs):
            t_min = ceil(limit_exprs[depth][0]())
            t_max = floor(limit_exprs[depth][1]())
            if t_max < t_min:
                return

            for ti in (t_min, t_max):
                t[depth] = ti
                param = param_exprs[depth]()
                recurse(depth+1, f'{desc} {t_min}≤{ti}≤{t_max},{letters[depth]}={param}')
        else:
            param = param_exprs[depth]()
            print(f'{desc} {letters[depth]}={param}')

    recurse()

This produces

 0≤0≤150,a=0 0≤0≤150,b=0 0≤150≤150,c=150 -150≤-150≤0,d=0 0≤0≤150,e=0 30≤30≤89,f=30 v=2999
 0≤0≤150,a=0 0≤0≤150,b=0 0≤150≤150,c=150 -150≤-150≤0,d=0 0≤0≤150,e=0 30≤89≤89,f=89 v=49
 0≤0≤150,a=0 0≤0≤150,b=0 0≤150≤150,c=150 -150≤-150≤0,d=0 0≤150≤150,e=150 -150≤-150≤-136,f=0 v=749
 0≤0≤150,a=0 0≤0≤150,b=0 0≤150≤150,c=150 -150≤-150≤0,d=0 0≤150≤150,e=150 -150≤-136≤-136,f=14 v=49
 0≤150≤150,a=150 -150≤-150≤0,b=0 0≤0≤150,c=0 0≤0≤150,d=0 0≤0≤150,e=0 66≤66≤125,f=66 v=2999
 0≤150≤150,a=150 -150≤-150≤0,b=0 0≤0≤150,c=0 0≤0≤150,d=0 0≤0≤150,e=0 66≤125≤125,f=125 v=49
 0≤150≤150,a=150 -150≤-150≤0,b=0 0≤0≤150,c=0 0≤0≤150,d=0 0≤150≤150,e=150 -150≤-150≤-100,f=0 v=2549
 0≤150≤150,a=150 -150≤-150≤0,b=0 0≤0≤150,c=0 0≤0≤150,d=0 0≤150≤150,e=150 -150≤-100≤-100,f=50 v=49
 0≤150≤150,a=150 -150≤-150≤0,b=0 0≤150≤150,c=150 -150≤-150≤0,d=0 0≤150≤150,e=150 -69≤-69≤-10,f=81 v=2999
 0≤150≤150,a=150 -150≤-150≤0,b=0 0≤150≤150,c=150 -150≤-150≤0,d=0 0≤150≤150,e=150 -69≤-10≤-10,f=140 v=49

An exhaustive generation is really not practical. At the very first viable value of the first parameter, and only the sixth viable value of the second parameter, the file takes up 500+MB of space and 26,000,000+ lines of text; so a full generation will easily consume "orders of magnitude above your disc size". You would be vastly better off making an "explorer interface" that understands the parameter bounds, and is able to navigate through slices of the parameter space based on user input. You have a seven-dimensional polytope; you need to make some concessions about how it's represented.

Beyond Brute Force

An 'explorer interface' could look like the following:

  • do not bother solving the Diophantine problem directly
  • allow the user to choose successively more fixed values for your parameters
  • do so by using multiple runs of a MIP (mixed-integer linear programming) library such as the venerable glpk

A demonstration:

import swiglpk as lp

from typing import Iterator


LETTERS = 'abcdefv'
N = len(LETTERS)


def check_lp(problem, code: int):
    msg = ''

    generic, primal, dual = (
        getattr(lp, f'glp_get_{fun}')(problem)
        for fun in ('status', 'prim_stat', 'dual_stat')
    )

    if generic != lp.GLP_OPT or primal != lp.GLP_FEAS or dual != lp.GLP_FEAS:
        statuses = {getattr(lp, f'GLP_{stat}'): stat
                    for stat in ('OPT', 'FEAS', 'INFEAS', 'NOFEAS', 'UNBND', 'UNDEF')}
        msg += ('Mehrotra feasibility constraints failed: '
                f'generic {statuses[generic]}, '
                f'primal {statuses[primal]}, '
                f'dual {statuses[dual]}')

    if code != 0:
        codes = {getattr(lp, k): k
                 for k in dir(lp) if k.startswith('GLP_E')}
        msg += f'glpk returned {codes[code]}'

    if msg:
        raise ValueError(msg)


def get_bounds(fixed: dict[str, int], letter: str) -> Iterator[float]:
    problem = lp.glp_create_prob()
    try:
        lp.glp_set_prob_name(problem, f'bounds_{letter}')
        lp.glp_set_obj_name(problem, 'single_parameter_range')
        '''
        1 auxiliary variable: the value of 
             -42a + 75b - 30c + 80d + 25e + 50f + v, which must equal -1
        7 structural variables: the values of "letters"
        '''
        lp.glp_add_rows(problem, 1)
        lp.glp_add_cols(problem, N)

        v_idx = lp.intArray(N+1)
        v_coeffs = lp.doubleArray(N+1)

        for j, (lett, v_coeff) in enumerate(zip(LETTERS, (-42, 75, -30, 80, 25, 50, 1)), 1):
            lp.glp_set_col_name(problem, j, lett)
            lp.glp_set_col_kind(problem, j, lp.GLP_IV)
            c = 1 if lett == letter else 0
            lp.glp_set_obj_coef(problem, j, c)

            if lett in fixed:
                lp.glp_set_col_bnds(problem, j, lp.GLP_FX, fixed[lett], fixed[lett])
            else:
                lp.glp_set_col_bnds(problem, j, lp.GLP_DB, 0, 3000 if lett=='v' else 150)

            v_idx[j] = j
            v_coeffs[j] = v_coeff

        lp.glp_set_mat_row(problem, 1, N, v_idx, v_coeffs)
        lp.glp_set_row_name(problem, 1, 'v_constraint')
        lp.glp_set_row_bnds(problem, 1, lp.GLP_FX, -1, -1)

        for direction in (lp.GLP_MIN, lp.GLP_MAX):
            lp.glp_set_obj_dir(problem, direction)

            level = lp.GLP_MSG_ERR  # lp.GLP_MSG_ON

            parm = lp.glp_smcp()
            lp.glp_init_smcp(parm)
            parm.msg_lev = level
            check_lp(problem, lp.glp_simplex(problem, parm))

            parm = lp.glp_iocp()
            lp.glp_init_iocp(parm)
            parm.msg_lev = level
            check_lp(problem, lp.glp_intopt(problem, parm))

            value = lp.glp_mip_col_val(problem, 1+LETTERS.index(letter))
            yield value
    finally:
        lp.glp_delete_prob(problem)


def main():
    fixed: dict[str, int] = {}

    while True:
        try:
            for letter in LETTERS:
                fixed_at = fixed.get(letter)
                if fixed_at is None:
                    bounds = get_bounds(fixed, letter)
                    print(f'{next(bounds):>3.0f} ≤ {letter}', end=' ', flush=True)
                    print(f'≤ {next(bounds):.0f}', flush=True)
                else:
                    print(f'{fixed_at:>3.0f} = {letter}')
        except ValueError as e:
            print(str(e))
        print()

        choice = input(f'Enter one of {LETTERS} to fix, '
                       f'any of {LETTERS.upper()} to free, or q to quit: ')
        if choice.lower() == 'q':
            break
        if choice.isupper():
            for c in choice:
                fixed.pop(c.lower(), None)
        else:
            index = LETTERS.find(choice.lower())
            if index >= 0:
                fixed[choice] = int(input(f'Enter new value of {choice}: '))
            print('Fixed parameters:', ', '.join(fixed.keys()))


if __name__ == '__main__':
    main()
  0 ≤ a ≤ 150
  0 ≤ b ≤ 143
  0 ≤ c ≤ 150
  0 ≤ d ≤ 134
  0 ≤ e ≤ 150
  0 ≤ f ≤ 150
  0 ≤ v ≤ 3000

Enter one of abcdefv to fix, any of ABCDEFV to free, or q to quit: d
Enter new value of d: 140
Mehrotra feasibility constraints failed: generic NOFEAS, primal NOFEAS, dual INFEAS

Enter one of abcdefv to fix, any of ABCDEFV to free, or q to quit: D
  0 ≤ a ≤ 150
  0 ≤ b ≤ 143
  0 ≤ c ≤ 150
  0 ≤ d ≤ 134
  0 ≤ e ≤ 150
  0 ≤ f ≤ 150
  0 ≤ v ≤ 3000

Enter one of abcdefv to fix, any of ABCDEFV to free, or q to quit: d
Enter new value of d: 130
141 ≤ a ≤ 150
  0 ≤ b ≤ 5
137 ≤ c ≤ 150
130 = d
  0 ≤ e ≤ 15
  0 ≤ f ≤ 7
  0 ≤ v ≤ 399

Enter one of abcdefv to fix, any of ABCDEFV to free, or q to quit: v
Enter new value of v: 200
148 ≤ a ≤ 148
  0 ≤ b ≤ 0
147 ≤ c ≤ 147
130 = d
  1 ≤ e ≤ 1
  0 ≤ f ≤ 0
200 = v
\$\endgroup\$
5
  • \$\begingroup\$ Thanks for going to such length to answer my question, but as you've mentioned writing a proper solver from scratch is neither easy nor desirable. That said, the nature of my question is not how to make a solver, but how to find the best solutions for a set of equations within certain boundaries. Despite being not scalable, brute-force approach is much easier to make and all I needed is a little bit of performance optimization to be able to run my program once and record the answers I need. After that I will have a ready to use table and the rest of the program will be of no use to me. \$\endgroup\$
    – Noby
    Jan 15, 2023 at 23:50
  • \$\begingroup\$ Ok, but... you shouldn't store a table. You should store the bounds of the solution space parameters. \$\endgroup\$
    – Reinderien
    Jan 15, 2023 at 23:51
  • \$\begingroup\$ I was planning on transfering it into an excell sheet later and just using raw values and some kind of simple algorithm to find the solutions. I'm not exactly sure how I would do that using your method \$\endgroup\$
    – Noby
    Jan 15, 2023 at 23:53
  • 1
    \$\begingroup\$ I can edit my answer to demonstrate. But this calls for the classic question: do you really need every answer? "What are you really doing?" \$\endgroup\$
    – Reinderien
    Jan 15, 2023 at 23:55
  • \$\begingroup\$ I'm just making a list of combinations for a specific interaction in a game. It's pretty useful to have such a list instead of re-calculating it every time \$\endgroup\$
    – Noby
    Jan 15, 2023 at 23:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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