yetanothermathprogrammingconsultant.blogspot.com Open in urlscan Pro
2a00:1450:4001:80e::2001  Public Scan

URL: http://yetanothermathprogrammingconsultant.blogspot.com/
Submission: On May 06 via api from GB — Scanned from GB

Form analysis 2 forms found in the DOM

Name: contact-form

<form name="contact-form">
  <p></p> Name <br>
  <input class="contact-form-name" id="ContactForm1_contact-form-name" name="name" size="30" type="text" value="">
  <p></p> Email <span style="font-weight: bolder;">*</span>
  <br>
  <input class="contact-form-email" id="ContactForm1_contact-form-email" name="email" size="30" type="text" value="">
  <p></p> Message <span style="font-weight: bolder;">*</span>
  <br>
  <textarea class="contact-form-email-message" cols="25" id="ContactForm1_contact-form-email-message" name="email-message" rows="5"></textarea>
  <p></p>
  <input class="contact-form-button contact-form-button-submit" id="ContactForm1_contact-form-submit" type="button" value="Send">
  <p></p>
  <div style="text-align: center; max-width: 222px; width: 100%">
    <p class="contact-form-error-message" id="ContactForm1_contact-form-error-message"></p>
    <p class="contact-form-success-message" id="ContactForm1_contact-form-success-message"></p>
  </div>
</form>

http://yetanothermathprogrammingconsultant.blogspot.com/search

<form action="http://yetanothermathprogrammingconsultant.blogspot.com/search" class="gsc-search-box" target="_top">
  <table cellpadding="0" cellspacing="0" class="gsc-search-box">
    <tbody>
      <tr>
        <td class="gsc-input">
          <input autocomplete="off" class="gsc-input" name="q" size="10" title="search" type="text" value="">
        </td>
        <td class="gsc-search-button">
          <input class="gsc-search-button" title="search" type="submit" value="Search">
        </td>
      </tr>
    </tbody>
  </table>
</form>

Text Content

YET ANOTHER MATH PROGRAMMING CONSULTANT

I am a full-time consultant and provide services related to the design,
implementation and deployment of mathematical programming, optimization and
data-science applications. I also teach courses and workshops. Usually I cannot
blog about projects I am doing, but there are many technical notes I'd like to
share. Not in the least so I have an easy way to search and find them again
myself. You can reach me at erwin@amsterdamoptimization.com.





WEDNESDAY, MAY 4, 2022


EVIL SUDOKU



On sudoku.com[1] there is now an "evil" difficulty level. I found these very
difficult to solve by hand (I usually give up).








 



Of course, a good MIP presolver/preprocessor devours these in no time. It should
solve models like this in less than 0.1 seconds using zero Simplex iterations
and zero B&B nodes.




Version identifier: 20.1.0.1 | 2021-04-07 | 3a818710c
CPXPARAM_Advance                                 0
CPXPARAM_Threads                                 1
CPXPARAM_MIP_Display                             4
CPXPARAM_MIP_Pool_Capacity                       0
CPXPARAM_MIP_Tolerances_AbsMIPGap                0
Generic callback                                 0x50
Tried aggregator 2 times.
MIP Presolve eliminated 168 rows and 580 columns.
Aggregator did 28 substitutions.
Reduced MIP has 128 rows, 122 columns, and 501 nonzeros.
Reduced MIP has 122 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.06 sec. (0.93 ticks)
Found incumbent of value 0.000000 after 0.06 sec. (1.66 ticks)

Root node processing (before b&c):
  Real time             =    0.06 sec. (1.67 ticks)
Sequential b&c:
  Real time             =    0.00 sec. (0.00 ticks)
                          ------------
Total (root+branch&cut) =    0.06 sec. (1.67 ticks)

--- MIP status (101): integer optimal solution.
--- Cplex Time: 0.06sec (det. 1.67 ticks)

Proven optimal solution
MIP Solution:            0.000000    (0 iterations, 0 nodes)
Final Solve:             0.000000    (0 iterations)

Best possible:           0.000000
Absolute gap:            0.000000
Relative gap:            0.000000





CBC is a bit terser (not always the case):




COIN-OR CBC      38.2.1 96226ea8 Feb 19, 2022          WEI x86 64bit/MS Window

COIN-OR Branch and Cut (CBC Library 2.10)
written by J. Forrest

Calling CBC main solution routine...
No integer variables - nothing to do

Solved to optimality (within gap tolerances optca and optcr).
MIP solution:    0.000000e+00   (0 nodes, 0.019 seconds)

Best possible:   0.000000e+00
Absolute gap:    0.000000e+00   (absolute tolerance optca: 0)
Relative gap:    0.000000e+00   (relative tolerance optcr: 0.0001)








REFERENCES






 1. https://sudoku.com/evil/




APPENDIX: GAMS MODEL




$ontext

   Evil Sudoku from sudoku.com

   References:
  
https://yetanothermathprogrammingconsultant.blogspot.com/2022/05/evil-sudoku.html

$offtext

set
  a      'all areas (rows,columns,squares)' /r1*r9,c1*c9,s1*s9/
  i(a)   'rows' /r1*r9/
  j(a)   'columns' /c1*c9/
  s(a)   'squares' /s1*s9/
  u(a,i,j) 'areas with unique values a=area name, (i,j) are the cells'
  k      'values' /1*9/;
;

*
* populate u(a,i,j)
*
u(i,i,j) = yes;
u(j,i,j) = yes;
u(s,i,j)$(ord(s)=3*ceil(ord(i)/3)+ceil(ord(j)/3)-3) = yes;
display u;

*
* given values
*
table v0(i,j)
   c1 c2 c3 c4 c5 c6 c7 c8 c9
r1                 9     5
r2  7                    1
r3        8  2        7     9
r4  6              4  9     8
r5     4     1
r6                          5
r7        2     3
r8  8              1  6     4
r9                       7
;

*
* MIP model
*
binary variable x(i,j,k);
x.fx(i,j,k)$(v0(i,j)=k.val) = 1;
variable z 'dummy objective';

equations
   dummy 'objective'
   unique(a,k) 'all-different'
   values(i,j) 'one value per cell'
;

dummy.. z =e= 0;
unique(a,k).. sum(u(a,i,j), x(i,j,k)) =e= 1;
values(i,j).. sum(k, x(i,j,k)) =e= 1;

model sudoku /all/;

*
* solve
*
solve sudoku minimizing z using mip;

*
* display solution
*
parameter v(i,j) 'solution';
v(i,j) = sum(k, k.val*x.l(i,j,k));
option v:0;
display v;







Most solvers have constraints for rows, columns, and squares. In this
formulation, I generalize this to "areas". The multidimensional set
u(a,i,j)u(a,i,j) contains for each area aa the corresponding cells (i,j)(i,j). 
You can display this set to view it. Once we have the set properly populated,
the uniqueness constraints become exceedingly simple. The idea is that
constraints are more difficult to debug than sets/parameters. So we spend more
effort in creating somewhat complex sets, but the pay-off is that the
constraints become very simple.


As usual, we have xi,j,k={1if cell (i,j) has value k0otherwisexi,j,k={1if
cell (i,j) has value k0otherwise To recover the Sodoku solution, we form
vi,j=∑kk⋅x∗i,j,kvi,j=∑kk⋅xi,j,k∗ 


The final solution is:





----     73 PARAMETER v  solution

            c1          c2          c3          c4          c5          c6          c7          c8          c9

r1           3           2           4           7           1           9           8           5           6
r2           7           5           9           8           6           3           4           1           2
r3           1           6           8           2           4           5           7           3           9
r4           6           7           1           3           5           4           9           2           8
r5           9           4           5           1           8           2           3           6           7
r6           2           8           3           9           7           6           1           4           5
r7           4           9           2           6           3           7           5           8           1
r8           8           3           7           5           2           1           6           9           4
r9           5           1           6           4           9           8           2           7           3



Posted by Erwin Kalvelagen at 4:21 PM 3 comments:





SUNDAY, MAY 1, 2022


GETTING RID OF NON-LINEARITIES


In [1], an attempt was made to implement the following non-linear constraint:
∑tmax{∑iCFi,t⋅qi−Lt,0}∑tLt≤0.015∑tmax{∑iCFi,t⋅qi−Lt,0}∑tLt≤0.015 Obviously code
like:

 

cfm = quicksum( max(quicksum(cf[i][t] * q[i] - L[t] for i in range(I)), 0) for t in range(T) /  quicksum(L[t] for t in range(T)) <= 0.015




is really not what we want to see. It is best, not to directly start up your
editor and write Python code like this. Rather, get a piece of paper, and focus
on the math.


There are two non-linearities: a division and a max()max() function. Note that I
am not sure which symbols are variables or parameters.  I'll try to explain some
possible cases below.





DIVISION



Exogenous divisions of the form xaxa where aa is a constant, are not an issue.
This is just really a linear coefficient 1/a1/a. If you have x∑iaix∑iai in the
constraints, I would advice to calculate α:=∑iaiα:=∑iai in advance, so we can
simplify things to xαxα It is a good habit to keep constraints as simple as
possible, as they are more difficult to debug than parameter assignments.


If we have something like xyxy (endogenous division) we need to be careful and
worry about division by zero (or, more generally, that the expression can assume
very large values if y≈0y≈0). If yy is a positive variable, we can protect
things a bit by using a lower bound, say y≥0.0001y≥0.0001. This is difficult to
do if yy is a free variable (we want a hole in the middle). If we have 
x∑iyix∑iyi it may be better to write this as: xzz=∑iyiz≥0.0001xzz=∑iyiz≥0.0001
at the expense of an extra variable and constraint. It has the additional
advantage of making the model less non-linear (fewer non-linear variables,
smaller non-linear part of the Jacobian). Note that we assumed here that
∑iyi≥0∑iyi≥0. If this expression can assume both negative and positive values,
this will not not work.


If the division is not embedded in other non-linear expressions, we often can
multiply things out. E.g. if we have the constraint:  x∑iyi≤0.015x∑iyi≤0.015
then we can write: x≤0.015⋅∑iyix≤0.015⋅∑iyi This is now linear. This is what can
apply to our original
constraint:∑tmax{∑iCFi,t⋅qi−Lt,0}≤0.015⋅∑tLt∑tmax{∑iCFi,t⋅qi−Lt,0}≤0.015⋅∑tLt





MAX FUNCTION



We can feed max{x,0}max{x,0} to a general-purpose NLP solver. However, this is
dangerous. Although this function is continuous, it is non-differentiable at
=0=0.


In the special case max{x,0}≤ymax{x,0}≤y we can write: x≤y0≤yx≤y0≤y If yy
represents an expression, you may want to use an extra variable to prevent
duplication.  


If we have ∑imax{xi,0}≤y∑imax{xi,0}≤y we can do: zi≥0zi≥xi∑izi≤yzi≥0zi≥xi∑izi≤y
Note that the variables zizi are a bit difficult to interpret. They can be
larger than max{xi,0}max{xi,0} if the inequality is not tight. The reason I use
zi≥max{xi,0}zi≥max{xi,0} is that it is sufficient for our purpose and that
zi=max{xi,0}zi=max{xi,0} is just way more difficult to represent.


This is the reformulation we can apply to our original problem:
zt≥∑iCFi,t⋅qi−Ltzt≥0∑tzt≤0.015⋅∑tLtzt≥∑iCFi,t⋅qi−Ltzt≥0∑tzt≤0.015⋅∑tLt This is
now completely linear: we have removed the division and the max(.,0)max(.,0)
function.





REFERENCES



 1. Piecewise Linear Function as
    Constraint, https://stackoverflow.com/questions/72026948/piecewise-linear-function-as-constraint




Posted by Erwin Kalvelagen at 1:48 PM No comments:





THURSDAY, APRIL 28, 2022


INVERTING A LARGE, DENSE MATRIX


In this post I show some experiences with inverting a large, dense matrix. The
variability in performance was a bit surprising to me. 



BACKGROUND



I have available in a GAMS GDX file a large matrix AA. This is an input-output
matrix. The goal is to find the so-called Leontief inverse (I−A)−1(I−A)−1 and
return this in a GAMS GDX for subsequent use [1]. It is well-known that instead
of calculating the inverse it is better to form an LU factorization of I−AI−A
and then use that for different rhs vectors bb (essentially a pair of easy
triangular solves). [2] However, economists are very much used to having access
to an explicit inverse.




A GAMS representation of the problem can look like:



set
  r 'regions' /r1*r2/
  p 'products' /p1*p3/
  rp(r,p)
;

rp(r,p) = yes;

table A(r,p,r,p)


        r1.p1  r1.p2  r1.p3  r2.p1  r2.p2  r2.p3
r1.p1    0.8    0.1
r1.p2    0.1    0.8    0.1
r1.p3           0.1    0.8    0.1
r2.p1                  0.1    0.8    0.1
r2.p2                                0.8    0.2
r2.p3                  0.1    0.1    0.1    0.7

;

parameters
   I(r,p,r,p)   'Identity matrix'
   IA(r,p,r,p)   'I-A'
;
alias(rp,rp2,rp3);
I(rp,rp) = 1;
IA(rp,rp2) = I(rp,rp2) - A(rp,rp2);
option I:1:2:2,IA:1:2:2;
display I,IA;


*-------------------------------------------------------------------
* inverse in GAMS
*-------------------------------------------------------------------

variable
   invIA(r,p,r,p)   'inv(I-A) calculated by a CNS model'
;

equation LinEq(r,p,r,p)   'IA*invIA=I';

LinEq(rp,rp2)..  sum(rp3, IA(rp,rp3)*invIA(rp3,rp2)) =e= I(rp,rp2);

model Linear /linEq/;
solve linear using cns;

option invIA:3:2:2;
display invIA.l;

* accuracy check.
* same test as in tool
scalar err 'max error in | (I-A)*inv(I-A) - I |';
err = smax((rp,rp2),abs(sum(rp3, IA(rp,rp3)*invIA.l(rp3,rp2)) - I(rp,rp2)));
display err;



The problem at hand has a 4-dimensional matrix AA (with both regions and sectors
as rows and columns). Of course, the real problem is very large. Because of
this, a pure GAMS model is not really feasible. In our real data set we have 306
regions and 63 sectors, which gives us a 19,278×19,27819,278×19,278 matrix. The
above model would have a variable and a constraint for each element in this
matrix. We would face a model with 371 million variables and constraints. Below,
we try to calculate the inverse using a standalone C++ program that reads a GDX
file with the AA matrix (and underlying sets for regions and sectors) and writes
a new GDX file with the Leontief inverse. 




ATTEMPT 1: MICROSOFT VISUAL C++ WITH THE EIGEN LIBRARY.



Here I use the Micosoft Visual C++ compiler, with the Eigen [3] library. Eigen
is a C++ template library. It is very easy to use (it is a header only library).
And it provides some interesting features:
 * It contains quite a few linear algebra algorithms, including LU
   factorization, solving, and inverting matrices.
 * It allows to write easy-to-read matrix expressions. E.g. to check if things
   are ok, I added a simple-minded accuracy check:
   max|(I−A)−1⋅(I−A)−I|max|(I−A)−1⋅(I−A)−I|  should be small. Here |.||.|
   indicates elementwise absolute values. This can be written as:
      double maxvalue = (INV * IA - MatrixXd::Identity(matrixSize,
   matrixSize)).array().abs().maxCoeff();
 * Support for parallel computing using OpenMP [4]. MSVC++ also supports this.



The results are:




Notes:
 * By default OMP (and Eigen) uses 16 threads. This machine has 8 cores and 16
   hardware threads.
 * Reading the input GDX file is quite fast: 1.3 seconds. (Note: this is read in
   "raw" mode for performance reasons). There is also some hashing going on.
 * The nonzero count of I−AI−A is 4.9 million. The inverse will be dense, so
   192782=372192782=372 million elements. We use here dense storage schemes and
   dense algorithms.
 * Instead of calling inverse(), I split the calculation explicitly into two
   parts: LU decomposition and solve. The time to perform the LU factorization
   looks reasonable to me, but the solve time is rather large. 
 * We can write the GDX file a bit faster by using compression. The writing time
   goes down from 29.4 to 14.8 seconds.



Why is the solve so slow? The following pictures may show the reason:





CPU usage during LU decomposition










CPU usage during solve





Basically, during the solve phase, we only keep about 2 cores at work. I would
expect that multiple rhs vectors (and a lot of them!) would give more
opportunity for parallelization. I have seen remarks that this operation is
memory bound, so adding more threads does not really help.



ATTEMPT 2: INTEL/CLANG WITH THE MKL LIBRARY.



Here I changed my toolset to Intel's Clang compiler icx, and its optimized MKL
library. The Clang compiler is part of the LLVM project, and is known for
excellent code generation. The MKL library is a highly optimized and tuned
version of the well-known LAPACK and BLAS libraries. It is possible to use MKL
as a back-end for Eigen, but here I just coded against MKL directly (I just have
one LU decomposition, an invert, and a matrix-multiplication in the check). So
how does this look like?





CPU Usage during solve






This looks much better. 


Notes:
 * MKL only uses 8 threads by default: the number of physical cores.
 * Both the LU decomposition, but especially the solve part is much faster.
 * Here we used compression when writing the GDX file.
 * I used Row Major Order memory layout (i.e. C style row-wise storage).
 * MKL is a bit more low-level. Basically, we use BLAS[6] and LAPACK[7].




CONCLUSION



Intel/Clang/MKL can be much faster than MSVC/Eigen on some linear algebra
operations. The total time to invert this matrix went from > 1000 seconds to  <
100 seconds:



Inverting a 19278x19278 matrix LibraryLU factorizationSolve MSVC++/Eigen49
seconds1047 seconds Intel/Clang/MKL24 seconds68 seconds






REFERENCES



 1. Matrix operations via GAMS Embedded Python: multi-regional Input-Output
    tables, https://yetanothermathprogrammingconsultant.blogspot.com/2021/08/matrix-operations-via-gams-embedded.html
 2. John D. Cook, Don't invert that
    matrix, https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/
 3. Eigen, https://eigen.tuxfamily.org/index.php
 4. OpenMP, https://www.openmp.org/.
 5. The LLVM Compiler Infrastructure, https://llvm.org/
 6. Basic Linear Algebra
    Subprograms, https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms.
 7. LAPACK, https://en.wikipedia.org/wiki/LAPACK.

 

Posted by Erwin Kalvelagen at 6:32 AM No comments:





SATURDAY, APRIL 16, 2022


EXPANDING VISUAL STUDIO IN THE TASK MANAGER



 






This is scary. Note that is in an idle state (nothing compiling, running, or
debugging at the moment). Soon we all need 128 GB of RAM.

Posted by Erwin Kalvelagen at 6:03 PM 1 comment:





THURSDAY, APRIL 14, 2022


GAMS: UNDOCUMENTED PUT FORMATS



The GAMS documentation mentions f.nr=1 (standard notation) and
f.nr=2 (scientific notation) for formatting numeric items using the PUT
statement. There are however a few more. Here we print values using different
formats:




       f.nr=1       f.nr=2       f.nr=3       f.nr=4
         1.20     1.20E+00 1.2000000000          1.2
         1.23     1.23E+00 1.2345000000       1.2345
 123450000.00     1.23E+08 123450000.00    123450000
         0.00     1.20E-07 0.0000001200       1.2E-7
         0.00     1.23E-07 0.0000001235    1.2345E-7






It looks like formats 3 and 4 are inspired by SAS PUT formats f12.10 and
best12.10.



REFERENCES

 1. Dictionary of
    formats, https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/leforinforref/p0z62k899n6a7wn1r5in6q5253v1.htm




APPENDIX: GAMS CODE




set
   i   'values'   /i1*i5/
   fmt 'formats'  /'f.nr=1'*'f.nr=4'/
;

parameter p(i) 'test values' /
   i1  1.2
   i2  1.2345
   i3  1.2345e8
   i4  1.2e-7
   i5  1.2345e-7
/;

*
* write all values in all formats
*

file f /put.txt/;
* right align text
f.lj = 1;
put f;
* write headers
loop(fmt,
   put fmt.tl:13
);
put /;
* write values
loop(i,
   loop(fmt,
      put ' ':0;
      f.nr=ord(fmt);
      put p(i);
   );
   put /;
);




Posted by Erwin Kalvelagen at 2:41 AM No comments:





SUNDAY, APRIL 10, 2022


REWRITING OLD GAMS CODE



It is always a good idea to revisit existing GAMS code and see if we can improve
it. Here is an example of an actual model.

The problem is that we want to set up a mapping set between two sets based on
the first two characters. If they are the same, the pair should be added to the
mapping. The old code looked like:




sets
   r 'regions' /
       AP   'Appalachia'
       CB   'Corn Belt'
       LA   'Lake States'
       /
   s 'subregions' /
       APN0501
       APN0502
       APN0504
       CBL0704
       CBM0404
       LAM0404
       LAM0406
       /
   rs(r,s) 'mapping regions<->subregions'
;


* old code:
PARAMETER value;
LOOP(R,  value('1',R) = 100*ord( R.tl,1) + ord( R.tl,2); );
LOOP(S,  value('2',S) = 100*ord( S.tl,1) + ord( S.tl,2); );
RS(R,S)$(value('1',R) = value('2',S)) = YES;


display value;
display rs;

  

To verify this indeed generates the correct result for set rs, we look at the
output:




----     27 PARAMETER value  

           AP          CB          LA     APN0501     APN0502     APN0504     CBL0704     CBM0404     LAM0404

1    6580.000    6766.000    7665.000
2                                        6580.000    6580.000    6580.000    6766.000    6766.000    7665.000

+     LAM0406

2    7665.000


----     28 SET rs  mapping regions<->subregions

       APN0501     APN0502     APN0504     CBL0704     CBM0404     LAM0404     LAM0406

AP         YES         YES         YES
CB                                             YES         YES
LA                                                                     YES         YES




The set rs is obviously correctly populated. The auxiliary parameter value
contains the ASCII values: e.g. AP has value 6580 as the decimal ASCII codes for
A and P are 65 and 80.


Notes:
 * The function ord(string, pos) returns the ASCII value of the character at
   position pos of the string.
 * The suffix .TL is the text string of the set element (internally GAMS does
   not work with the strings but rather an integer id for each set element; here
   we make it explicit we want the string).
 * The loops are not really needed. We could have used:
   value('1',R) = 100*ord( R.tl,1) + ord( R.tl,2);
   value('2',S) = 100*ord( S.tl,1) + ord( S.tl,2);
 * The mapping set is a very powerful concept in GAMS. It is somewhat like a
   dict in Python, except it is not one way. It can map from r to s but also
   from s to r.



We can simplify the code to just:



* new code:
* compare first two characters in r and s
rs(r,s) = ord(r.tl,1)=ord(s.tl,1) and ord(r.tl,2)=ord(s.tl,2);
display rs;


This gives the same result:




----     22 SET rs  mapping regions<->subregions

       APN0501     APN0502     APN0504     CBL0704     CBM0404     LAM0404     LAM0406

AP         YES         YES         YES
CB                                             YES         YES
LA                                                                     YES         YES






Notes:
 * No need for the auxiliary parameter value.
 * No loops.
 * We could add a check that we have no unmapped subregions. E.g. by:
   abort$(card(rs)<>card(s)) "There are unmapped subregions";
 * This will also catch cases where the casing does not match.
 * Critical note: GAMS has very poor support for strings. This should have been
   rectified a long time ago. There is no excuse for this deficiency. 




CONCLUSION



It is always worthwhile to revisit existing models and clean them up. During
development, especially when under time pressure, we often are happy once our
code works. As the lifetime of models can be surprisingly long, and thus new
generations of users and modelers start working with old models, it is a good
idea to do a cleanup round to make the code as readable as possible.  

Posted by Erwin Kalvelagen at 7:20 PM No comments:





SATURDAY, MARCH 5, 2022


BOSTOCK'S COMMAND-LINE CARTOGRAPHY TUTORIAL UNDER WINDOWS



In [1], an excellent tutorial is presented about the process of making maps. It
is a little bit dated, so here I develop some Windows-based scripts that make it
possible to follow these tutorials. The goal is two-fold:



 1. Make things work for Windows. I am comfortable with the Unix command line,
    but, by far, most of my colleagues are not. To allow for easier sharing, I
    used Windows CMD.
 2. Update the commands so they all function again. Some URLs have changed a
    little bit, and we need to use a specific version of d3. If you want to use
    the original Unix commands, and you encounter some issues, you may want to
    check how I repaired them. 



The output format of these scripts is svg. SVG files are plain text, represent
vectors (instead of bitmaps), and can be read directly by a standard web
browser.

 

PART 1




The first thing to do is to install node.js if you don't have this available
already. Using [2] this process is quite straightforward. Our first batch file
will follow the steps in [1].




 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

::---------------------------------------------------------------------------
::
::   Mike Bostock
::   Command-Line Cartography, Part 1
::   A tour of d3-geo’s new command-line interface.
:: 
::   https://medium.com/@mbostock/command-line-cartography-part-1-897aa8f8ca2c
::
::   Needs node.js. Install from: https://nodejs.org/en/download/
::
::   we translate UNIX command line to Windows CMD
::---------------------------------------------------------------------------

:: download shape file from Census Bureau (region 06 = California)
::
:: curl 'http://www2.census.gov/geo/tiger/GENZ2014/shp/cb_2014_06_tract_500k.zip' -o cb_2014_06_tract_500k.zip
:: use https instead, curl is part of windows
curl -O https://www2.census.gov/geo/tiger/GENZ2014/shp/cb_2014_06_tract_500k.zip 

:: unzip -o cb_2014_06_tract_500k.zip
:: use tar instead (part of windows)
tar -xf cb_2014_06_tract_500k.zip

:: install node.js package
:: and extract geojson
call npm install -g shapefile
call shp2json cb_2014_06_tract_500k.shp -o ca.json

:: perform projection 
call npm install -g d3-geo-projection
call geoproject "d3.geoConicEqualArea().parallels([34, 40.5]).rotate([120, 0]).fitSize([960, 960], d)" < ca.json > ca-albers.json

:: create svg
call geo2svg -w 960 -h 960 < ca-albers.json > ca-albers.svg

:: launch browser
start ca-albers.svg





Notes:



 * The curl command is available on newer versions of Windows.
 * Windows does not have unzip, but it has tar which can unzip files.
 * Use https instead of HTTP.
 * The use of quotes is different under Windows compared to Unix. 



The result should be the following vector image displayed in your browser:






Read more »
Posted by Erwin Kalvelagen at 2:29 AM No comments:





FRIDAY, FEBRUARY 18, 2022


WANTED: BETTER ERROR MESSAGES


In [1] a user tries the following constraint:

> mdl.addConstrs((t[i,k] * X[i,j,k] - te1[i] <= 5) >> (z1[i,k] == 1) for i,j,k
> in arcos if i != 0 and i != 23)

(t and X and variables). This is not accepted, and the following error message
is issued:



> AttributeError: 'gurobipy.QuadExpr' object has no attribute 'getVar'



This is a terrible error message. It does not describe the real, underlying
problem, and it does not tell the user what to do next. Developers typically do
not pay much attention to error messages. Often they seem more aimed to make
life easy for the developer than providing assistance to an already confused
user. In this case, I believe the modeling tool did not even issue an error
message, but rather made assumptions without checking, and left it to the Python
run-time system to complain. 


This constraint has two problems: 
 * Indicator constraints have a simple binary variable as condition, not a
   general (boolean) expression
 * Only linear constraints are supported in indicator constraints



Maybe a better message is:



> Indicator constraints have the form:
> 
> 
>         binary variable == 0 ==> linear constraint
> 
> 
> or
> 
> 
>         binary variable == 1 ==> linear constraint 
> 
> 
> Your constraint does not have this form. You need to reformulate it to fit
> this format, or use a different (big-M) formulation.




Well-formed and meaningful error messages are very important. The user is
already confused, otherwise (s)he would formulate a correct constraint. Adding
confusing error messages to this confusion is not a good idea. A good error
message could have prevented this post. We can think of error messages as the UI
(User Interface) when it really matters: when the user needs help. Developers
should pay much more attention to this.

A question for another day is of course: why are only linear constraints
supported?




REFERENCES






 1. AttributeError: 'gurobipy.QuadExpr' object has no attribute
    'getVar', https://stackoverflow.com/questions/71153625/attributeerror-gurobipy-quadexpr-object-has-no-attribute-getvar/71158372








Posted by Erwin Kalvelagen at 11:06 AM 4 comments:





WEDNESDAY, FEBRUARY 16, 2022


SLSQP ORIGINAL PAPER



This is somewhat hard to find (although I did), so I share it. This is the
original paper[1] describing the SLSQP solver used in scipy[2]:


Read more »
Posted by Erwin Kalvelagen at 8:00 AM No comments:





TUESDAY, FEBRUARY 15, 2022


VISUALIZATION: ANIMATING FLOW ALONG AN EDGE



In a large network, it is a bit of a pain to visualize the flow. Here is an
attempt to show the shortest path in a random sparse directed graph. The
shortest path is of course a min-cost flow of one unit. I used cytoscape.js [1]
to generate the picture.















Read more »
Posted by Erwin Kalvelagen at 10:29 AM No comments:



Older Posts Home

Subscribe to: Posts (Atom)



ABOUT ME

Erwin Kalvelagen View my complete profile



HOME PAGE

Amsterdam Optimization Modeling Group LLC




CONTACT FORM



Name




Email *




Message *











BLOG ARCHIVE

 * ▼  2022 (15)
   * ▼  May (2)
     * Evil Sudoku
     * Getting rid of non-linearities
   * ►  April (4)
   * ►  March (1)
   * ►  February (4)
   * ►  January (4)

 * ►  2021 (52)
   * ►  December (2)
   * ►  November (3)
   * ►  October (3)
   * ►  September (4)
   * ►  August (3)
   * ►  July (5)
   * ►  June (5)
   * ►  May (2)
   * ►  April (4)
   * ►  March (9)
   * ►  February (6)
   * ►  January (6)

 * ►  2020 (60)
   * ►  December (7)
   * ►  November (5)
   * ►  October (5)
   * ►  September (4)
   * ►  August (5)
   * ►  July (4)
   * ►  June (6)
   * ►  May (7)
   * ►  April (4)
   * ►  March (5)
   * ►  February (2)
   * ►  January (6)

 * ►  2019 (50)
   * ►  December (5)
   * ►  November (5)
   * ►  October (5)
   * ►  September (2)
   * ►  August (5)
   * ►  July (3)
   * ►  June (4)
   * ►  May (5)
   * ►  April (4)
   * ►  March (2)
   * ►  February (4)
   * ►  January (6)

 * ►  2018 (79)
   * ►  December (4)
   * ►  November (4)
   * ►  October (6)
   * ►  September (4)
   * ►  August (5)
   * ►  July (7)
   * ►  June (4)
   * ►  May (11)
   * ►  April (11)
   * ►  March (9)
   * ►  February (6)
   * ►  January (8)

 * ►  2017 (94)
   * ►  December (8)
   * ►  November (6)
   * ►  October (8)
   * ►  September (8)
   * ►  August (3)
   * ►  July (6)
   * ►  June (5)
   * ►  May (11)
   * ►  April (11)
   * ►  March (12)
   * ►  February (7)
   * ►  January (9)

 * ►  2016 (166)
   * ►  December (15)
   * ►  November (12)
   * ►  October (18)
   * ►  September (11)
   * ►  August (12)
   * ►  July (6)
   * ►  June (13)
   * ►  May (21)
   * ►  April (13)
   * ►  March (13)
   * ►  February (17)
   * ►  January (15)

 * ►  2015 (80)
   * ►  December (15)
   * ►  November (8)
   * ►  October (7)
   * ►  September (1)
   * ►  August (3)
   * ►  July (4)
   * ►  June (8)
   * ►  May (3)
   * ►  April (8)
   * ►  March (8)
   * ►  February (6)
   * ►  January (9)

 * ►  2014 (45)
   * ►  December (6)
   * ►  November (6)
   * ►  October (2)
   * ►  September (6)
   * ►  August (2)
   * ►  July (3)
   * ►  June (3)
   * ►  May (1)
   * ►  April (7)
   * ►  March (2)
   * ►  February (5)
   * ►  January (2)

 * ►  2013 (72)
   * ►  December (1)
   * ►  November (1)
   * ►  October (3)
   * ►  September (2)
   * ►  August (4)
   * ►  July (6)
   * ►  June (9)
   * ►  May (10)
   * ►  April (5)
   * ►  March (7)
   * ►  February (9)
   * ►  January (15)

 * ►  2012 (86)
   * ►  December (12)
   * ►  November (7)
   * ►  October (10)
   * ►  September (4)
   * ►  August (10)
   * ►  July (7)
   * ►  June (2)
   * ►  May (4)
   * ►  April (7)
   * ►  March (11)
   * ►  February (7)
   * ►  January (5)

 * ►  2011 (76)
   * ►  December (9)
   * ►  November (10)
   * ►  October (10)
   * ►  September (7)
   * ►  August (4)
   * ►  July (9)
   * ►  June (6)
   * ►  May (6)
   * ►  April (6)
   * ►  March (1)
   * ►  February (6)
   * ►  January (2)

 * ►  2010 (89)
   * ►  December (5)
   * ►  November (2)
   * ►  October (10)
   * ►  September (4)
   * ►  August (7)
   * ►  July (8)
   * ►  June (10)
   * ►  May (9)
   * ►  April (8)
   * ►  March (10)
   * ►  February (7)
   * ►  January (9)

 * ►  2009 (194)
   * ►  December (12)
   * ►  November (10)
   * ►  October (11)
   * ►  September (17)
   * ►  August (9)
   * ►  July (15)
   * ►  June (18)
   * ►  May (22)
   * ►  April (21)
   * ►  March (29)
   * ►  February (14)
   * ►  January (16)

 * ►  2008 (197)
   * ►  December (16)
   * ►  November (9)
   * ►  October (24)
   * ►  September (16)
   * ►  August (16)
   * ►  July (30)
   * ►  June (48)
   * ►  May (38)




SEARCH THIS BLOG







Awesome Inc. theme. Powered by Blogger.



Diese Website verwendet Cookies von Google, um Dienste anzubieten und Zugriffe
zu analysieren. Deine IP-Adresse und dein User-Agent werden zusammen mit
Messwerten zur Leistung und Sicherheit für Google freigegeben. So können
Nutzungsstatistiken generiert, Missbrauchsfälle erkannt und behoben und die
Qualität des Dienstes gewährleistet werden.Weitere InformationenOk