Skip to content

Commit fbe19ea

Browse files
kabdelhakkabdelhak
authored
[NB] implement sparsity coloring (#8711)
- implement greedy distance 2 algorithm from 'what color is your jacobian?' publication - currently only uses column compression - idea for row compression also implemented -> currently not used but computed and printed out (constant time so not a deal) -> does not consider indirect row dependencies through inner equations, needs to be added!
1 parent 2dbc7f7 commit fbe19ea

File tree

7 files changed

+190
-95
lines changed

7 files changed

+190
-95
lines changed

‎OMCompiler/Compiler/NBackEnd/Classes/NBackendDAE.mo‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@ public
4242
import Events = NBEvents;
4343
import NFFlatten.FunctionTree;
4444
import Jacobian = NBJacobian;
45-
import NBJacobian.SparsityPattern;
45+
import NBJacobian.{SparsityPattern, SparsityColoring};
4646
import StrongComponent = NBStrongComponent;
4747
import NBSystem;
4848
import NBSystem.System;
@@ -154,7 +154,7 @@ public
154154
for i in 1:arrayLength(bdae.comps) loop
155155
tmp := tmp + StrongComponent.toString(bdae.comps[i], i) + "\n";
156156
end for;
157-
tmp := tmp + SparsityPattern.toString(bdae.sparsityPattern, bdae.sparsityColoring);
157+
tmp := tmp + SparsityPattern.toString(bdae.sparsityPattern) + "\n" + SparsityColoring.toString(bdae.sparsityColoring);
158158
then tmp;
159159

160160
case HESSIAN() then StringUtil.headline_1("Hessian: " + str) + "\n" +

‎OMCompiler/Compiler/NBackEnd/Modules/3_Post/NBJacobian.mo‎

Lines changed: 137 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -175,7 +175,7 @@ public
175175
VarData varData;
176176
EqData eqData;
177177
SparsityPattern sparsityPattern;
178-
list<list<ComponentRef>> sparsityColoring = {};
178+
SparsityColoring sparsityColoring = SparsityColoring.lazy(EMPTY_SPARSITY_PATTERN, UnorderedMap.new<CrefLst>(ComponentRef.hash, ComponentRef.isEqual));
179179
algorithm
180180

181181
if listLength(jacobians) == 1 then
@@ -187,8 +187,6 @@ public
187187
local
188188
VarData tmpVarData;
189189
SparsityPattern tmpPattern;
190-
Integer size1, size2;
191-
list<list<ComponentRef>> coloring1, coloring2;
192190

193191
case BackendDAE.JACOBIAN(varData = tmpVarData as VarData.VAR_DATA_JAC(), sparsityPattern = tmpPattern) algorithm
194192
jacType := jac.jacType;
@@ -210,20 +208,7 @@ public
210208
seed_vars := listAppend(tmpPattern.seed_vars, seed_vars);
211209
partial_vars := listAppend(tmpPattern.partial_vars, partial_vars);
212210
nnz := nnz + tmpPattern.nnz;
213-
214-
// combine the sparsity colorings since all are independent
215-
size1 := listLength(sparsityColoring);
216-
size2 := listLength(jac.sparsityColoring);
217-
(coloring1, coloring2) := if size1 > size2
218-
then (sparsityColoring, jac.sparsityColoring)
219-
else (jac.sparsityColoring, sparsityColoring);
220-
221-
// fill up the smaller coloring with empty groups
222-
for i in 1:intAbs(size1-size2) loop
223-
coloring2 := {} :: coloring2;
224-
end for;
225-
226-
sparsityColoring := List.threadMap(coloring1, coloring2, listAppend);
211+
sparsityColoring := SparsityColoring.combine(sparsityColoring, jac.sparsityColoring);
227212
then ();
228213

229214
else algorithm
@@ -297,8 +282,11 @@ public
297282
end match;
298283
end jacobianTypeString;
299284

300-
type SparsityPatternCol = tuple<ComponentRef, list<ComponentRef>> "residual, {independents}";
301-
type SparsityPatternRow = SparsityPatternCol "independent, {residuals}";
285+
// necessary as wrapping value type for UnorderedMap
286+
type CrefLst = list<ComponentRef>;
287+
288+
type SparsityPatternCol = tuple<ComponentRef, CrefLst> "partial_vars, {seed_vars}";
289+
type SparsityPatternRow = SparsityPatternCol "seed_vars, {partial_vars}";
302290

303291
uniontype SparsityPattern
304292
record SPARSITY_PATTERN
@@ -311,7 +299,6 @@ public
311299

312300
function toString
313301
input SparsityPattern pattern;
314-
input SparsityColoring coloring;
315302
output String str = StringUtil.headline_2("Sparsity Pattern (nnz: " + intString(pattern.nnz) + ")");
316303
protected
317304
ComponentRef cref;
@@ -333,42 +320,29 @@ public
333320
str := str + "(" + ComponentRef.toString(cref) + ")\t depends on:\t" + ComponentRef.listToString(dependencies) + "\n";
334321
end for;
335322
end if;
336-
if colEmpty and rowEmpty then
337-
str := str + "\n<empty sparsity pattern>\n";
338-
end if;
339-
str := str + "\n" + StringUtil.headline_3("Sparsity Coloring Groups");
340-
if not listEmpty(coloring) then
341-
for group in coloring loop
342-
str := str + ComponentRef.listToString(group) + "\n";
343-
end for;
344-
else
345-
str := str + "<empty sparsity coloring>\n";
346-
end if;
347323
end toString;
348324

349-
// necessary as wrapping value type for UnorderedMap
350-
type CrefLst = list<ComponentRef>;
351-
352325
function create
353326
input VariablePointers seedCandidates;
354327
input VariablePointers partialCandidates;
355328
input Option<array<StrongComponent>> strongComponents "Strong Components";
356329
input JacobianType jacType;
357330
output SparsityPattern sparsityPattern;
358331
output SparsityColoring sparsityColoring;
332+
protected
333+
UnorderedMap<ComponentRef, CrefLst> map;
359334
algorithm
360-
sparsityPattern := match strongComponents
335+
(sparsityPattern, map) := match strongComponents
361336
local
362337
array<StrongComponent> comps;
363338
list<ComponentRef> seed_vars, seed_vars_array, partial_vars, tmp;
364-
UnorderedMap<ComponentRef, list<ComponentRef>> map;
365339
UnorderedSet<ComponentRef> set;
366340
list<SparsityPatternCol> cols = {};
367341
list<SparsityPatternRow> rows = {};
368342
Integer nnz = 0;
369343

370344
case SOME(comps) guard(arrayEmpty(comps)) algorithm
371-
then EMPTY_SPARSITY_PATTERN;
345+
then (EMPTY_SPARSITY_PATTERN, UnorderedMap.new<CrefLst>(ComponentRef.hash, ComponentRef.isEqual));
372346

373347
case SOME(comps) algorithm
374348
// get all relevant crefs
@@ -420,7 +394,7 @@ public
420394
(_, tmp) := col;
421395
nnz := nnz + listLength(tmp);
422396
end for;
423-
then SPARSITY_PATTERN(cols, rows, seed_vars, partial_vars, nnz);
397+
then (SPARSITY_PATTERN(cols, rows, seed_vars, partial_vars, nnz), map);
424398

425399
case NONE() algorithm
426400
Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because of missing strong components."});
@@ -433,30 +407,144 @@ public
433407
end match;
434408

435409
// create coloring
436-
sparsityColoring := createEmptyColoring(sparsityPattern);
410+
sparsityColoring := SparsityColoring.PartialD2ColoringAlg(sparsityPattern, map);
437411

438412
if Flags.isSet(Flags.DUMP_SPARSE) then
439-
print(toString(sparsityPattern, sparsityColoring));
413+
print(toString(sparsityPattern) + "\n" + SparsityColoring.toString(sparsityColoring));
440414
end if;
441415
end create;
442416

443417
function createEmpty
444418
output SparsityPattern sparsityPattern = EMPTY_SPARSITY_PATTERN;
445-
output SparsityColoring sparsityColoring = createEmptyColoring(sparsityPattern);
419+
output SparsityColoring sparsityColoring = EMPTY_SPARSITY_COLORING;
446420
end createEmpty;
421+
end SparsityPattern;
422+
423+
constant SparsityPattern EMPTY_SPARSITY_PATTERN = SPARSITY_PATTERN({}, {}, {}, {}, 0);
424+
constant SparsityColoring EMPTY_SPARSITY_COLORING = SPARSITY_COLORING(listArray({}), listArray({}));
425+
426+
type SparsityColoringCol = list<ComponentRef> "seed variable lists belonging to the same color";
427+
type SparsityColoringRow = SparsityColoringCol "partial variable lists for each color (multiples allowed!)";
428+
429+
uniontype SparsityColoring
430+
record SPARSITY_COLORING
431+
"column wise coloring with extra row sparsity information"
432+
array<SparsityColoringCol> cols;
433+
array<SparsityColoringRow> rows;
434+
end SPARSITY_COLORING;
435+
436+
function toString
437+
input SparsityColoring sparsityColoring;
438+
output String str = StringUtil.headline_2("Sparsity Coloring");
439+
protected
440+
Boolean empty = arrayLength(sparsityColoring.cols) == 0;
441+
algorithm
442+
if empty then
443+
str := str + "\n<empty sparsity pattern>\n";
444+
end if;
445+
for i in 1:arrayLength(sparsityColoring.cols) loop
446+
str := str + "Color (" + intString(i) + ")\n"
447+
+ " - Column: " + ComponentRef.listToString(sparsityColoring.cols[i]) + "\n"
448+
+ " - Row: " + ComponentRef.listToString(sparsityColoring.rows[i]) + "\n\n";
449+
end for;
450+
end toString;
447451

448-
function createEmptyColoring
449-
"creates an empty coloring that just groups each independent variable individually"
452+
function lazy
453+
"creates a lazy coloring that just groups each independent variable individually
454+
and implies dependence for each row"
450455
input SparsityPattern sparsityPattern;
451-
output SparsityColoring sparsityColoring = {};
456+
input UnorderedMap<ComponentRef, CrefLst> map;
457+
output SparsityColoring sparsityColoring;
458+
protected
459+
array<SparsityColoringCol> cols;
460+
array<SparsityColoringRow> rows;
461+
algorithm
462+
cols := listArray(list({cref} for cref in sparsityPattern.seed_vars));
463+
rows := arrayCreate(arrayLength(cols), sparsityPattern.partial_vars);
464+
sparsityColoring := SPARSITY_COLORING(cols, rows);
465+
end lazy;
466+
467+
function PartialD2ColoringAlg
468+
"author: kabdelhak 2022-03
469+
taken from: 'What Color Is Your Jacobian? Graph Coloring for Computing Derivatives'
470+
https://doi.org/10.1137/S0036144504444711
471+
A greedy partial distance-2 coloring algorithm. Slightly adapted to also track row sparsity."
472+
input SparsityPattern sparsityPattern;
473+
input UnorderedMap<ComponentRef, CrefLst> map;
474+
output SparsityColoring sparsityColoring;
475+
protected
476+
array<ComponentRef> cref_lookup;
477+
UnorderedMap<ComponentRef, Integer> index_lookup;
478+
array<Boolean> color_exists;
479+
array<Integer> coloring, forbidden_colors;
480+
array<list<ComponentRef>> col_coloring, row_coloring;
481+
Integer color;
482+
list<SparsityColoringCol> cols_lst = {};
483+
list<SparsityColoringRow> rows_lst = {};
452484
algorithm
453-
sparsityColoring := list({cref} for cref in sparsityPattern.seed_vars);
454-
end createEmptyColoring;
485+
// integer to cref and reverse lookup arrays
486+
cref_lookup := listArray(sparsityPattern.seed_vars);
487+
index_lookup := UnorderedMap.new<Integer>(ComponentRef.hash, ComponentRef.isEqual, Util.nextPrime(listLength(sparsityPattern.seed_vars)));
488+
for i in 1:arrayLength(cref_lookup) loop
489+
UnorderedMap.add(cref_lookup[i], i, index_lookup);
490+
end for;
455491

456-
end SparsityPattern;
492+
// create empty colorings
493+
coloring := arrayCreate(arrayLength(cref_lookup), 0);
494+
forbidden_colors := arrayCreate(arrayLength(cref_lookup), 0);
495+
color_exists := arrayCreate(arrayLength(cref_lookup), false);
496+
col_coloring := arrayCreate(arrayLength(cref_lookup), {});
497+
row_coloring := arrayCreate(arrayLength(cref_lookup), {});
498+
499+
for i in 1:arrayLength(cref_lookup) loop
500+
for row_var /* w */ in UnorderedMap.getSafe(cref_lookup[i], map) loop
501+
for col_var /* x */ in UnorderedMap.getSafe(row_var, map) loop
502+
color := coloring[UnorderedMap.getSafe(col_var, index_lookup)];
503+
if color > 0 then
504+
forbidden_colors[color] := i;
505+
end if;
506+
end for;
507+
end for;
508+
(_, color) := Array.findFirstOnTrueWithIdx(forbidden_colors, function intNe(i2 = i));
509+
coloring[i] := color;
510+
// also save all row dependencies of this color
511+
row_coloring[color] := listAppend(row_coloring[color], UnorderedMap.getSafe(cref_lookup[i], map));
512+
color_exists[color] := true;
513+
end for;
457514

458-
constant SparsityPattern EMPTY_SPARSITY_PATTERN = SPARSITY_PATTERN({}, {}, {}, {}, 0);
459-
type SparsityColoring = list<list<ComponentRef>> "list of independent variable groups belonging to the same color";
515+
for i in 1:arrayLength(coloring) loop
516+
col_coloring[coloring[i]] := cref_lookup[i] :: col_coloring[coloring[i]];
517+
end for;
518+
519+
// traverse in reverse to have correct ordering in the end)
520+
for i in arrayLength(color_exists):-1:1 loop
521+
if color_exists[i] then
522+
cols_lst := col_coloring[i] :: cols_lst;
523+
rows_lst := row_coloring[i] :: rows_lst;
524+
end if;
525+
end for;
526+
527+
sparsityColoring := SPARSITY_COLORING(listArray(cols_lst), listArray(rows_lst));
528+
end PartialD2ColoringAlg;
529+
530+
function combine
531+
"combines sparsity patterns by just appending them because they are supposed to
532+
be entirely independent of each other."
533+
input SparsityColoring coloring1;
534+
input SparsityColoring coloring2;
535+
output SparsityColoring coloring_out;
536+
protected
537+
SparsityColoring smaller_coloring;
538+
algorithm
539+
// append the smaller to the bigger
540+
(coloring_out, smaller_coloring) := if arrayLength(coloring2.cols) > arrayLength(coloring1.cols) then (coloring2, coloring1) else (coloring1, coloring2);
541+
542+
for i in 1:arrayLength(smaller_coloring.cols) loop
543+
coloring_out.cols[i] := listAppend(coloring_out.cols[i], smaller_coloring.cols[i]);
544+
coloring_out.rows[i] := listAppend(coloring_out.rows[i], smaller_coloring.rows[i]);
545+
end for;
546+
end combine;
547+
end SparsityColoring;
460548

461549
protected
462550
// ToDo: all the DAEMode stuff is probably incorrect!

‎OMCompiler/Compiler/NSimCode/NSimJacobian.mo‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -360,7 +360,7 @@ public
360360
protected
361361
list<Integer> tmp;
362362
algorithm
363-
for group in listReverse(coloring) loop
363+
for group in listReverse(arrayList(coloring.cols)) loop
364364
try
365365
tmp := list(SimVar.getIndex(UnorderedMap.getSafe(cref, sim_map)) for cref in group);
366366
else

‎testsuite/simulation/modelica/NBackend/array_handling/simple_der_for.mos‎

Lines changed: 7 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -135,19 +135,13 @@ simulate(simple_der_for); getErrorString();
135135
// ($DER.x[2]) depends on: {x[2]}
136136
// ($DER.x[1]) depends on: {x[1]}
137137
//
138-
// Sparsity Coloring Groups
139-
// ----------------------------------------
140-
// {x[1]}
141-
// {x[2]}
142-
// {x[3]}
143-
// {x[4]}
144-
// {x[5]}
145-
// {x[6]}
146-
// {x[7]}
147-
// {x[8]}
148-
// {x[9]}
149-
// {x[10]}
150-
// {x[11]}
138+
// ========================================
139+
// Sparsity Coloring
140+
// ========================================
141+
// Color (1)
142+
// - Column: {x[11],x[10],x[9],x[8],x[7],x[6],x[5],x[4],x[3],x[2],x[1]}
143+
// - Row: {$DER.x[1],$DER.x[2],$DER.x[4],$DER.x[5],$DER.x[6],$DER.x[7],$DER.x[8],$DER.x[9],$DER.x[10],$DER.x[11]}
144+
//
151145
//
152146
// record SimulationResult
153147
// resultFile = "simple_der_for_res.mat",

‎testsuite/simulation/modelica/NBackend/functions/builtin_functions.mos‎

Lines changed: 7 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -321,22 +321,13 @@ simulate(builtin_functions); getErrorString();
321321
// ($DER.x[2]) depends on: {}
322322
// ($DER.x[1]) depends on: {}
323323
//
324-
// Sparsity Coloring Groups
325-
// ----------------------------------------
326-
// {x[1]}
327-
// {x[2]}
328-
// {x[3]}
329-
// {x[4]}
330-
// {x[5]}
331-
// {x[6]}
332-
// {x[7]}
333-
// {x[8]}
334-
// {x[9]}
335-
// {x[10]}
336-
// {x[11]}
337-
// {x[12]}
338-
// {x[13]}
339-
// {x[14]}
324+
// ========================================
325+
// Sparsity Coloring
326+
// ========================================
327+
// Color (1)
328+
// - Column: {x[14],x[13],x[12],x[11],x[10],x[9],x[8],x[7],x[6],x[5],x[4],x[3],x[2],x[1]}
329+
// - Row: {$DER.x[14]}
330+
//
340331
//
341332
// record SimulationResult
342333
// resultFile = "builtin_functions_res.mat",

‎testsuite/simulation/modelica/NBackend/minimalTearing.mos‎

Lines changed: 14 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -74,9 +74,13 @@ simulate(minimalTearing); getErrorString();
7474
// ----------------------------------------
7575
// ($RES_SIM_0) depends on: {y}
7676
//
77-
// Sparsity Coloring Groups
78-
// ----------------------------------------
79-
// {y}
77+
// ========================================
78+
// Sparsity Coloring
79+
// ========================================
80+
// Color (1)
81+
// - Column: {y}
82+
// - Row: {$RES_SIM_0}
83+
//
8084
//
8185
// ################################################################################
8286
//
@@ -130,9 +134,13 @@ simulate(minimalTearing); getErrorString();
130134
// ----------------------------------------
131135
// ($RES_SIM_0) depends on: {y}
132136
//
133-
// Sparsity Coloring Groups
134-
// ----------------------------------------
135-
// {y}
137+
// ========================================
138+
// Sparsity Coloring
139+
// ========================================
140+
// Color (1)
141+
// - Column: {y}
142+
// - Row: {$RES_SIM_0}
143+
//
136144
//
137145
// record SimulationResult
138146
// resultFile = "minimalTearing_res.mat",

0 commit comments

Comments
 (0)