Implementing simplex algorithm in lua

[surreal epoch]

@gleaming crane Made a thread 😄

[gleaming crane]

Do you have a single file for the functions, or a zip?

[surreal epoch]

So I'm trying to implement the simplex algorithm in lua. I'm following the algorithm via a description on wikipedia: https://en.wikipedia.org/wiki/Simplex_algorithm#Simplex_tableau

[surreal epoch]

It's messy because I was in the middle of debugging with a bunch of logs

[gleaming crane]

If you use VS Code, you can add breakpoints

[surreal epoch]

I was doing that but the debugger has been running too slowly for my purposes

Right now I could probably add it back in, but at one point I was building a dependency graph on all recipes in pyanodons and randomizing it while keeping the game finishable and well, that took a while even without the debugger going

[gleaming crane]

That's a lot

[surreal epoch]

Yeah it's fun stuff 😄

I heard py devs don't use the debugger either

[gleaming crane]

The only thing I can think right now is a number of locals for these tables could make it easier to read

[surreal epoch]

For performance reasons as well

You mean performance wise?

[gleaming crane]

I mean easier to read

[surreal epoch]

Oh duh

On the topic of being bad at reading...

You mean instead of doing params.whatever just do whatever

[gleaming crane]

One sec

[surreal epoch]

Something that complicates things is that I'll really know it works when it starts working on the recipes/items in vanilla, but my interface between vanilla stuff and the corresponding linear program could also be wrong

Right now it's saying that it takes zero recipe crafts to start from ore and make chemical science pack, so I know something is going awry

[gleaming crane]

    local m = params.row
    local n = params.col
    local matrix = params.matrix

    local scale_factor = 1 / matrix[m][n]
    for j=1,#matrix[m] do
        matrix[m][j] = scale_factor * matrix[m][j]
    end

That isn't the best example, but it makes it easier for *me* to follow

[surreal epoch]

nods

Would it help if I cleaned up the code with some comments too?

Might take a min

[gleaming crane]

I'm picky about comments. I try to be a little strict with my style choices lol

[surreal epoch]

ah okay

[gleaming crane]

I was thinking 'comments are for 1) exceptions and 2) why something "isn't easy to follow" or "doesn't look like what how you think it should look" '

They're basically the same thing, if you have a "comment" it should tell me something really important, though stating steps, especially for an algorithm, makes sense

[surreal epoch]

Same, but thought I might leave some for what part of the wikipedia page it corresponds to

But at the same time I guess you can just cross-reference

[gleaming crane]

I won't tell you to not make comments, but make them useful

[surreal epoch]

No I'm aware of this philosophy

[gleaming crane]

One small suggestion, is after a line like:

    for i=2,#params.matrix do

since you're using a table, assign it to a local

    for i=2,#params.matrix do
        local row = params.matrix[i]

Right?

[surreal epoch]

Sure sure, I'll try that

[gleaming crane]

    local n = params.col

    local L = #params.matrix
    for i=2,L do
        local row = params.matrix[i]
        if row[n] >= 1 and row[L] / row[n] > max then
            max = row[L] / row[n]
            max_ind = i
        end
    end

now I can read it :p

[surreal epoch]

thanks so much for the help

I've been stewing over this for months, returning every now and then just to be confused

[gleaming crane]

I haven't helped with the testing, but this is making it easier to parse

[surreal epoch]

I guess it was a thanks in advance

[gleaming crane]

Is this function correct?

function find_entering_var (params)
    local ROWS = #params.matrix
    local COLUMNS = #params.matrix[1]
    for i=ROWS, COLUMNS-1 do
        local row = params.matrix[1]
        if row[i] > 0 then
            params.col = i
            return true
        end
    end

    return false
end

The numbers to use don't seem matched

#params.matrix is the number of rows, and #params.matrix[1] is the length of the first row (the number of columns)

And params.matrix[1] is always row one 🤔

function find_entering_var (params)
    local ROWS = #params.matrix
    local COLUMNS = #params.matrix[1]
    local row_one = params.matrix[1]
    for i=ROWS, COLUMNS-1 do
        if row_one[i] > 0 then
            params.col = i
            return true
        end
    end

    return false
end

[surreal epoch]

Yeah, I'm not sure how else to get the number of columns

The matrix will look like this, so I'm only interested in the (k+1)-th row, where k is the number of columns

[gleaming crane]

Getting the number of columns is OK, but this iteration doesn't make sense? Should it be going from 1-COLUMNS?

Ohhhhh I think I get it

[surreal epoch]

Though honestly I could probably iterate through the first few as well, it won't trigger the condition

[gleaming crane]

It only checks the upper triangular numbers

[surreal epoch]

Uh... something like that

It's not a square matrix

[gleaming crane]

Er, yeah, it doesn't do that right now

[surreal epoch]

So it checks the part after the square bit

[gleaming crane]

Oh ok ok I think I'm actually catching on

It's a smooshed together identity matrix and another one?

[surreal epoch]

yep

Yeah so the matrix should look like in the image each time, I want to find an element that's "good" and then I do a pivot operation with that element

[gleaming crane]

Interesting. Should it stop at COLUMNS or COLUMNS - 1? Lua is one indexed of course

[surreal epoch]

Basically, if the top row has a positive element ever (other than the last element), we are good to keep going

[surreal epoch]

Should be columns-1? I'll double check the algorithm, but like in the picture the last column is "special"

Good catches though

[gleaming crane]

Always comparing row one seems a little strange, but is it like, submatrices are taken?

I've seen some fancy matrix stuff in another mod (the Command and Conquer Tiberium one)

[surreal epoch]

I'm not sure if I'm inputting the problem to it correctly, but it keeps saying the optimal value is zero, even when it isn't hm

[gleaming crane]

My linear algebra knowledge is very limited

[surreal epoch]

It's actually straightforward, just the first row is special

[surreal epoch]

I don't know if it would help, it didn't really help me other than knowing what pivots were

According to wikipedia, we just stop when there are no positive values in the first row

[gleaming crane]

I see.

[surreal epoch]

Otherwise, we find a positive value and do something with it, that should be straightforward at least

That gives us a column, we find a row by following this rule that to me is black magic

[gleaming crane]

Do you have an example you have worked through by hand?

[surreal epoch]

Then we get an element of the matrix to do the pivot on

[surreal epoch]

I couldn't find a nontrivial problem online but I just did

[gleaming crane]

The easiest way is to compare your work vs what the code says each time

[surreal epoch]

So I'll try working through this

[gleaming crane]

Yeah

[surreal epoch]

If you look up linear programming problems you get two-variable deals...

And then the solution will just be like "graph it and look at the graph"

I found this though earlier and didn't realize it had solved examples: https://online-optimizer.appspot.com/?model=https:%2F%2Fraw.githubusercontent.com%2Fyozw%2Flio-files%2Fmaster%2Fchapter01%2Fdiet.mod

[gleaming crane]

I might just rewrite the first function a bit, so it's easier to follow? Or I could give you a bit and see what you find?

[surreal epoch]

Rewriting sounds good. The first function I have like 80% confidence in, which is pretty high compared to the others at least 😛

Do you understand what it does?

[gleaming crane]

I have an idea... I may see something

I think this rewrite is accurate

[surreal epoch]

Aha! the last row shouldn't be all zeroes

Something's wrong with this bit

[gleaming crane]

The first section could probably be rewritten to make more sense

[surreal epoch]

Okay, so this was added after there were problems but the fact that it's missing an absolute value probably doesn't help...

[gleaming crane]

If it's slightly positive or negative it is ok to round to zero, yes?

[surreal epoch]

Hopefully? I was having issues with numerical stability

Basically, if there's an entry that's 0.00000000001 and it stumbles on it later, it will turn into 10000000000 or smth I didn't count the zeroes

[gleaming crane]

Ah

If you want a name for that concept, I like to use "epsilon", or, "e"

[surreal epoch]

And the way the solutions are, at least in factorio, we shouldn't expect 10,000,000,000 to come up? I'm not too sure what the right magic number is there

Good wording

[gleaming crane]

ϵ

I don't use the symbol but yeah, in one of my mods I have a local e = 0.0001 or so at the top of a file

[surreal epoch]

That doesn't look right?

ε

Oh that's how it looks

Code epsilon is weird

[gleaming crane]

Weird, no, mine's slightly less rounded in the middle

... I don't know

[surreal epoch]

You're right

weird

[gleaming crane]

ε yours is right

I guess the other one is what it used in "TeX"?

or I don't know lol.

https://wumbo.net/symbols/epsilon/ this site goes back and forth

[surreal epoch]

Do you not know latex?

[gleaming crane]

ϵ Greek Lunate Epsilon Symbol[1],
ε Greek Small Letter Epsilon[1]

I have heard of it, but never used it. Don't know any of the jargon

[surreal epoch]

It's still saying the score is zero hrmm

But for the toy problem it correctly answered -126

However, the toy problem only ever needs to use the pivot function, it was a little too toy

FML the next example looks insanely complicated 🥺

[gleaming crane]

I did some more condensing, in solve_system, and near the end

That code does not look like Lua lol . Are some of those 'slices'?

[surreal epoch]

the "could"?

[gleaming crane]

Sorry, code

[surreal epoch]

What do you mean by slices?

And how doesn't it look like Lua 😂

[gleaming crane]

I guess they're not slices. Slices make subsets of the storage object in some languages

[surreal epoch]

That sentence went over my head

[gleaming crane]

Some people who come from math / other languages would ask like "how do I just make a new table from..." and it's like, "that's not something you really do in Lua..."

For a matrix, someone might want to make a new smaller matrix and operate on that, affecting the original matrix. Lua does not provide built-in ways to do that.

[surreal epoch]

Well, the only place I need an auxilliary matrix is when I'm swapping columns; I save one column, move the other on top of it, and then move the auxiliarry/copied column to the other's spot.

Whyyyyy😔

I might need to increase epsilon

[gleaming crane]

That's a negative number... how did it sneak through

[surreal epoch]

No no negative numbers are allowed

Just not close to zero

[gleaming crane]

Ah, that's why you were suggesting math.abs

[surreal epoch]

Ye

I didn't expect numerical stability to be such an issue lol

I might have to encode the entries as pairs of numbers with special addition/multiplication/etc. function so that I can do exact math...

[gleaming crane]

rational numbers?? What are we, Americans?

[surreal epoch]

I am 😂

[gleaming crane]

Same.

[surreal epoch]

Ah yes, my favorite number, "probably one half"

Although it's still saying the objective function is exactly zero, so I suspect this isn't (yet) the issue

NIL

[gleaming crane]

Heh

[surreal epoch]

My old enemy friend

[gleaming crane]

Oh no I made an error I think

No... huh that's actually what was written

[surreal epoch]

where?

[gleaming crane]

    if permutation == nil then
        permutation = {}
        for i=1, COLUMNS do
            permutation[i] = i
        end
        params.permutation = permutation
    end

I thought I missed an indexing operation in permutation[i] = i

Assigning some other matrix's [i]th element to permutation[i]

... well I did make a bigger error

anyway, you're working from your code

[surreal epoch]

ye

[surreal epoch]

oh haha

permutation is just a debugging thing for now but will be used later

It's so I can keep track of the column swaps

It's also how I figured out nothing was being swapped

[gleaming crane]

A small thing is I can't exactly rely on making some values local, particularly row and col, because they could change in a called function. Only you would know if they should be made local, or if locals need to be re-assigned after a function call.

Refactoring might solve any possible confusion (when using locals and parameters that might've changed), but I wouldn't think about that until I was too confused or it worked.

[surreal epoch]

row is used upstream of another function in that file even I believe

[gleaming crane]

Upvalues are dangerous, I'd never use them

[surreal epoch]

But yeah refactoring isn't as high on my priority list right now given that I at least understand everything

[surreal epoch]

Well, it's not an upvalue per se, it's included in the table that's being passed

Which is partly why I opted for params.row rather than just row

[gleaming crane]

Don't forgot I guess, Lua can return and assign multiple values

[surreal epoch]

I thought this was how you returned multiple values 😛

I remember I looked it up and it had to do with tables so I ended up doing this way since it seemed the same. Perhaps not.

Oh do you just separate via commas

[gleaming crane]

It's a difference in style, personally I'd keep values which can change separate from those which don't (clarified by arguments)

[surreal epoch]

return thing1, thing2, thing3

[gleaming crane]

Right

[surreal epoch]

Right I remember now

[gleaming crane]

If something like working row and working column were passed in as arguments, then the function return could include the new ones

[surreal epoch]

yep

By the way, I learned that the matrix I was passing in has more rows than columns... so it wasn't exactly feasible

Apparently there are more items in factorio than there are recipes

Wait have I shared the code that actually constructs the matrix to be passed yet?

Got rid of the all-zeroes rows (those without a recipe) and still not working

This is the code, it's in a larger (much larger) file that constructs a "dependency graph". The relevant bit is that dependency_graph[prg.get_key({name = [name of an item or a fluid], type = "itemorfluid"}] will get a table that has a "prereqs" member that's itself a table saying the recipes that make it and how much of it they make. The dependents member is the recipes that consume it.

Recipes are the same with item ingredients as prereqs and items produced as dependents.

The linear program is essentially "minimize the number of recipes run to produce 1 chemical science pack while making sure there are no negative amounts of items" and I include "recipes" for producing iron ore etc. out of thin air as well.

Wait a frickin second it did something

Does this look like the number of "recipes" you need for chemical science?..

Seems high to me

[gleaming crane]

uhhhhh

-5624? probably not

[surreal epoch]

Disregard the negative sign btw

Yeah I'd expect it to be in the realm of 300-700

No wait probably lower

[gleaming crane]

For unmodded, yeah

Even for crafting 3 packs, it'd be like, 50, 20, and 10 for engines, red circuits, and sulfur

Very loose estimate

[surreal epoch]

No wait I think it was numerical stability again

[gleaming crane]

A bit more if we count mining operations and plates which are 2x... uh oh

[surreal epoch]

So an automation science pack is 8

Right?

“Year” is “gear” sorry

[gleaming crane]

That makes sense

I can read cursive, fam

[surreal epoch]

well... it seems to work for automation science pakcs?

I guess oil processing is messing it up big time somehow?

[gleaming crane]

Are you ignoring barrels?

[surreal epoch]

This seems correct for logistic science packs

[surreal epoch]

Uhhhhh

Those don't sound like a great combination with numerical instability

Computer: "Wait, I can trade 1 empty barrel for 1.00000000002 empty barrels according to these calculations? Great!"

[surreal epoch]

Okay I did, but now the matrix doesn't have enough columns again (which corresponds to being overdetermined; there aren't enough recipes)

I added some void recipes to compensate but now... well... this is what logistic science packs look like:

I think it's due to numerical stability in large part... I'll have to fix that and I'm tired now 😦

[gleaming crane]

Same, it's very late

Well, early

200k for a logistic pack... I would've just ignored fill and empty barreling from the initial recipes, hope you didn't try anything more complicated for now

[surreal epoch]

Just gonna keep writing here with my progress... I realized even doing fractions isn't gonna extend to the future because recipes in general utilize probabilities