Remove global variables and use PrecompileTools #9

odow · 2025-04-22T23:53:45Z

I attempted to make a fairly minimal set of changes. ~~, but there's a difference between the two models I need to debug.~~

Before

% julia --project=. demo.jl 
Create demo: 21.260735 seconds (28.36 M allocations: 1.897 GiB, 3.32% gc time, 99.10% compilation time: <1% of which was recompilation)
Running HiGHS 1.10.0 (git hash: fd8665394e): Copyright (c) 2025 HiGHS under MIT licence terms
LP   has 604 rows; 950 cols; 1920 nonzeros
Coefficient ranges:
  Matrix [1e-02, 1e+04]
  Cost   [1e+00, 1e+00]
  Bound  [1e-01, 1e+01]
  RHS    [9e-02, 8e+03]
Presolving model
265 rows, 287 cols, 965 nonzeros  0s
Dependent equations search running on 25 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
228 rows, 245 cols, 826 nonzeros  0s
Presolve : Reductions: rows 228(-376); columns 245(-705); elements 826(-1094)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     1.2861637975e+04 Pr: 30(74.9799) 0s
        157     3.8706151237e+04 Pr: 0(0); Du: 0(4.9738e-14) 0s
Solving the original LP from the solution after postsolve
Model status        : Optimal
Simplex   iterations: 157
Objective value     :  3.8706151237e+04
Relative P-D gap    :  3.7595872395e-16
HiGHS run time      :          0.01
Solve model: 0.124531 seconds (57.21 k allocations: 4.288 MiB, 89.49% compilation time: 51% of which was recompilation)

This PR

(base) oscar@MacBookPro TIMES % julia --project=. demo.jl
Create demo: 1.204018 seconds (762.20 k allocations: 47.372 MiB, 1.88% gc time, 92.00% compilation time: 44% of which was recompilation)
Running HiGHS 1.10.0 (git hash: fd8665394e): Copyright (c) 2025 HiGHS under MIT licence terms
LP   has 604 rows; 950 cols; 1920 nonzeros
Coefficient ranges:
  Matrix [1e-02, 1e+04]
  Cost   [1e+00, 1e+00]
  Bound  [1e-01, 1e+01]
  RHS    [9e-02, 8e+03]
Presolving model
265 rows, 287 cols, 965 nonzeros  0s
Dependent equations search running on 25 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
228 rows, 245 cols, 826 nonzeros  0s
Presolve : Reductions: rows 228(-376); columns 245(-705); elements 826(-1094)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     1.2861637975e+04 Pr: 30(74.9799) 0s
        157     3.8706151237e+04 Pr: 0(0); Du: 0(4.9738e-14) 0s
Solving the original LP from the solution after postsolve
Model status        : Optimal
Simplex   iterations: 157
Objective value     :  3.8706151237e+04
Relative P-D gap    :  3.7595872395e-16
HiGHS run time      :          0.00
Solve model: 0.108321 seconds (57.21 k allocations: 4.288 MiB, 92.00% compilation time: 55% of which was recompilation)

odow · 2025-04-23T00:18:54Z

I think there's still 10-100X in performance lying around:

We could read in a proper data structure instead of the mix of data frames and dicts and named tuples and strings and symbols
The constraints could be tidied up. I'd love to see a link to the GAMS equivalent of this monster

JuMP.@constraint(
        model,
        EQG_COMBAL[(r, t, c, s) in indices["EQG_COMBAL"]],
        (
            !isnothing(data[:RHS_COMPRD]) && ((r, t, c, s) in data[:RHS_COMPRD]) ? ComPrd[(r, t, c, s)] :
            (
                sum(
                    (
                        (r, p, c) in data[:RPC_STG] ?
                        sum(
                            sum(
                                StgFlo[(r, v, t, p, c, ts, "OUT")] *
                                get(data[:RS_FR], (r, s, ts), 0) *
                                (
                                    1 + (
                                        !isnothing(data[:RTCS_FR]) ?
                                        get(data[:RTCS_FR], (r, t, c, s, ts), 0) : 0
                                    )
                                ) *
                                data[:STG_EFF][r, v, p] for v in get(sets.RTP_VNT, (r, t, p), Set())
                            ) for ts in sets.RPC_TS[r, p, c]
                        ) :
                        sum(
                            sum(
                                PrcFlo[(r, v, t, p, c, ts)] *
                                get(data[:RS_FR], (r, s, ts), 0) *
                                (
                                    1 + (
                                        !isnothing(data[:RTCS_FR]) ?
                                        get(data[:RTCS_FR], (r, t, c, s, ts), 0) : 0
                                    )
                                ) for v in get(sets.RTP_VNT, (r, t, p), Set());
                                init = 0,
                            ) for ts in sets.RPC_TS[r, p, c]
                        )
                    ) for p in get(sets.RCIO_P, (r, c, "OUT"), Set()) if (r, t, p) in data[:RTP_VARA];
                    init = 0,
                ) + sum(
                    sum(
                        sum(
                            IreFlo[(r, v, t, p, c, ts, "IMP")] *
                            get(data[:RS_FR], (r, s, ts), 0) *
                            (
                                1 + (
                                    !isnothing(data[:RTCS_FR]) ?
                                    get(data[:RTCS_FR], (r, t, c, s, ts), 0) : 0
                                )
                            ) for v in get(sets.RTP_VNT, (r, t, p), Set());
                            init = 0,
                        ) for ts in sets.RPC_TS[r, p, c]
                    ) for p in get(sets.RCIE_P, (r, c, "IMP"), Set()) if (r, t, p) in data[:RTP_VARA];
                    init = 0,
                )
            ) * data[:COM_IE][r, t, c, s]
        ) >=
        sum(
            (
                (r, p, c) in data[:RPC_STG] ?
                sum(
                    sum(
                        StgFlo[(r, v, t, p, c, ts, "IN")] *
                        get(data[:RS_FR], (r, s, ts), 0) *
                        (1 + (!isnothing(data[:RTCS_FR]) ? get(data[:RTCS_FR], (r, t, c, s, ts), 0) : 0))
                        for v in get(sets.RTP_VNT, (r, t, p), Set());
                        init = 0,
                    ) for ts in sets.RPC_TS[r, p, c]
                ) :
                (sum(
                    sum(
                        PrcFlo[(r, v, t, p, c, ts)] *
                        get(data[:RS_FR], (r, s, ts), 0) *
                        (1 + (!isnothing(data[:RTCS_FR]) ? get(data[:RTCS_FR], (r, t, c, s, ts), 0) : 0))
                        for v in get(sets.RTP_VNT, (r, t, p), Set());
                        init = 0,
                    ) for ts in sets.RPC_TS[r, p, c]
                ))
            ) for p in get(sets.RCIO_P, (r, c, "IN"), Set()) if (r, t, p) in data[:RTP_VARA];
            init = 0,
        ) +
        sum(
            sum(
                sum(
                    IreFlo[(r, v, t, p, c, ts, "EXP")] *
                    get(data[:RS_FR], (r, s, ts), 0) *
                    (1 + (!isnothing(data[:RTCS_FR]) ? get(data[:RTCS_FR], (r, t, c, s, ts), 0) : 0)) for
                    v in get(sets.RTP_VNT, (r, t, p), Set());
                    init = 0,
                ) for ts in sets.RPC_TS[r, p, c]
            ) for p in get(sets.RCIE_P, (r, c, "EXP"), Set()) if (r, t, p) in data[:RTP_VARA];
            init = 0,
        ) +
        get(data[:COM_PROJ], (r, t, c), 0) * data[:COM_FR][r, t, c, s]
    )

olejandro · 2025-04-23T10:12:22Z

Thanks a lot @odow! I have not gone through all the changes, but to answer your comments/questions briefly:

We could read in a proper data structure instead of the mix of data frames and dicts and named tuples and strings and symbols

What could that data structure be? Please bare in mind that the data that TIMES source code receives is different from the one used here, because the prototype does not include preprocessing. Or did you mean a data structure for the internal sets / parameters?

The constraints could be tidied up. I'd love to see a link to the GAMS equivalent of this monster

Here you go: https://github.com/etsap-TIMES/TIMES_model/blob/master/eqcombal.mod
There are probably some differences between the current version and the one that the prototype is based on.

olejandro · 2025-04-23T10:14:11Z

Project.toml


 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-GAMS = "1ca51c6a-1b4d-4546-9ae1-53e0a243ab12"


We should definitely keep this one to enable the users to use the solvers they may have through their GAMS license.

No, this should be removed from TIMES.jl. The user can separately do ] add GAMS and then set_optimizer(model, GAMS.Optimizer). There's no need to force every person to install GAMS.jl

I understand your point, but the rationale for including it would be similar to HiGHs. Many of the TIMES users have access to solvers through their perpetual GAMS licenses, so it should be easy for them to try this out.

src/compute_indexes.jl

olejandro · 2025-04-23T11:04:06Z

src/sets.jl

These ones should actually be sets (i.e unique elements, order doesn't matter, etc...). Why to change them to arrays?

We can revert if you like. I would have assumed that the input database did not accept duplicate values. If so, that seems like it's either a sanity check somewhere else, or an error here.

Also note that Set{String} is implemented as a Dict{String,Nothing}, so there is some cost as opposed to Vector{String}. It depends how large the list of elements is. Fewer than 100 and it's probably better to use Vector.

olejandro · 2025-04-23T11:35:27Z

src/load_data.jl

    if row_number > 0 && col_number == 1
        # One-dimensional set
-        value = Set(values(df[!, 1]))
+        return values(df[!, 1])


I believe this should be a set, i.e. elements should be unique.

olejandro · 2025-04-23T11:38:24Z

src/load_data.jl

+            return Dict(Tuple.(eachrow(df[:, DataFrames.Not(:value)])) .=> df.value)
        else
-            value = Set(Tuple.(eachrow(df)))
+            return Tuple.(eachrow(df))


Also here, the tuples should be unique.

olejandro · 2025-04-23T15:23:32Z

Interesting! When I run the code twice (by executing include('demo.jl')), the time it takes to create model on the second run with the changes seems the same / very similar to the original code (i.e. 0.3-0.4 seconds). Does this mean that the changes mostly improve the compilation time?

olejandro · 2025-04-23T16:31:12Z

src/TIMES.jl


+PrecompileTools.@setup_workload begin
+    PrecompileTools.@compile_workload begin
+        create_model("PROTO.db3")


Does this have a desired effect also when "PROTO.db3" changes?

It does't matter if the user's .db3 is different to PROTO.db3. All that matters is that PROTO.db3 is a small representative file that hits the correct code paths. You could make it even smaller.

Ok, understood. How about this one: https://julialang.github.io/PackageCompiler.jl?
Majority of the TIMES users do not interact with the code directly, but mostly focus on the input data (thanks to a very good separation of data and structure from the code). Would this be a good way to distribute the code for those users?

I wouldn't worry about how to distribute code yet. Let's assume you're >1 year away from delivering an end-user product. I assume the first users will be more technically savvy. They can manually install Julia.

For future, yes, you could do something like PackageCompiler, but there are number of gotchas. It wouldn't encourage you to spend time on this just yet. There's also a push within the core Julia developers to improve static compilation of binaries and apps like this. I would expect the landscape in a year to be different to today.

The goal of the project is to provide guidance to the community on the benefits and barriers to a possible migration to Julia/JuMP. To be useful, it should take into account how the TIMES source code is used and distributed at the moment. That is why it is an important question. It is good to know that things are being developed. 💪

There is a lot of work happening on juliac which aims at producing a binary from a Julia source code, but it's still in early stage of development. An essential trim feature will appear in Julia v1.12.

Thanks @blegat! It is really good to know where things are headed.

odow · 2025-04-24T00:30:18Z

Does this mean that the changes mostly improve the compilation time?

Yes. But it's also an artifact of the tiny problem size. For larger models, you'd probably see a difference between the two formulations.

odow · 2025-04-24T00:32:53Z

Here you go: https://github.com/etsap-TIMES/TIMES_model/blob/master/eqcombal.mod

That is actually horrific. I don't understand how it is possible to read or test that.

siddharth-krishna

Thanks Oscar, not a Julia programmer, but at a high level your changes look good to me. Very impressive speedup, and reduction in memory!

What data structure would you recommend for reading the input into?

src/compute_indexes.jl

odow · 2025-04-24T21:53:01Z

I would keep as much of the data munging in the DB as possible. Extract only the data that you need when you need it. Do all the filtering and joining when needed.

See, for example:

olejandro · 2025-04-25T08:39:03Z

Does this mean that the changes mostly improve the compilation time?

Yes. But it's also an artifact of the tiny problem size. For larger models, you'd probably see a difference between the two formulations.

:-) Will be aiming at getting a larger problem soon!

olejandro · 2025-04-25T08:40:32Z

I would keep as much of the data munging in the DB as possible. Extract only the data that you need when you need it. Do all the filtering and joining when needed.

See, for example:
* https://jump.dev/JuMP.jl/stable/tutorials/linear/multi/

* https://jump.dev/JuMP.jl/stable/tutorials/linear/multi_commodity_network/

Why so? Is it not much slower than in-memory operations?

olejandro · 2025-04-25T08:44:37Z

Here you go: https://github.com/etsap-TIMES/TIMES_model/blob/master/eqcombal.mod

That is actually horrific. I don't understand how it is possible to read or test that.

This version may be more comprehensible!

odow · 2025-04-27T20:44:51Z

Why so? Is it not much slower than in-memory operations?

Performance is not the main reason. You can write equally fast normal Julia code. Main reasons are:

To reduce the memory required to store various copies of the data
To have a single unified representation of the data, instead of dicts and sets and lists and tables and databases
Because you can probably rethink how to form the constraints instead of generating sparse sets of sets of sets and iterating over them. A common approach would be: do a groupby, add a new column to the table which represents the decision variable. Constraint is a sum product of two columns.

You may want to look at https://github.com/TulipaEnergy/TulipaEnergyModel.jl/tree/main/src for inspiration.

odow added 2 commits April 23, 2025 11:39

Rewrite prototype to remove global variables and use PrecompileTools

852b943

Fix

a85066a

olejandro reviewed Apr 23, 2025

View reviewed changes

src/compute_indexes.jl Show resolved Hide resolved

olejandro reviewed Apr 23, 2025

View reviewed changes

olejandro changed the title ~~Rewrite prototype to remove global variables and use PrecompileTools~~ Remove global variables and use PrecompileTools Apr 23, 2025

olejandro merged commit 7300b80 into olejandro:main Apr 23, 2025

olejandro reviewed Apr 23, 2025

View reviewed changes

siddharth-krishna reviewed Apr 24, 2025

View reviewed changes

src/compute_indexes.jl Show resolved Hide resolved

odow deleted the od/rewrite branch April 24, 2025 21:53

Remove global variables and use PrecompileTools #9

Remove global variables and use PrecompileTools #9

Uh oh!

Conversation

odow commented Apr 22, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Before

This PR

Uh oh!

odow commented Apr 23, 2025

Uh oh!

olejandro commented Apr 23, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

olejandro commented Apr 23, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

olejandro Apr 25, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

odow commented Apr 24, 2025

Uh oh!

odow commented Apr 24, 2025

Uh oh!

siddharth-krishna left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

odow commented Apr 24, 2025

Uh oh!

olejandro commented Apr 25, 2025

Uh oh!

olejandro commented Apr 25, 2025

Uh oh!

olejandro commented Apr 25, 2025

Uh oh!

odow commented Apr 27, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

odow commented Apr 22, 2025 •

edited

Loading

olejandro Apr 25, 2025 •

edited

Loading