Variable number of nested loops

roy437 · ‎12-16-2020

Hi,

How to write a Fortran code with a variable number of nested loops.

Thank you.

roy437 · ‎12-18-2020

I know the initial values of k() and l() that do not satisfy the given equation. I would like to find the solution of the equation that is closest to the initial solution (I would use the least squares method).

mecej4 · ‎12-18-2020

As you may have noticed, no objective function was specified in the model that I provided. You can try adding as objective for minimization the distance (euclidean? something else?) from the initial point of attraction.

roy437 · ‎12-18-2020

Very complicated for me, I don't know the program at all.

mecej4 · ‎12-18-2020

Here is the AMPL model file, augmented to minimize the Euclidean distance from an initial point k0, l0:

# Roy437 test problem, added requirement:
# {k, l} to be as close to {k0, l0} as possible

# Declarations
param M integer := 11;
param V integer := 268177985;
set IDX = 1..M;
var k {IDX} integer >= 0 <= 9;
var l {IDX} integer >= 0 <= 9;
param k0 {IDX} integer;
param l0 {IDX} integer;
param n {IDX} integer >= 1;
param j {IDX} integer >= 1;

# Objective function
minimize dist2:
   sum {i in IDX} ((k[i]-k0[i])^2 + (l[i]-l0[i])^2);

# Constraints
subj to C: sum{i in IDX} (100*j[i]+10*k[i]+l[i])*n[i] = V;

# Values of Parameters
data ;
param: n :=  1 811 2 263 3 715 4 3 5 534 6 286 7 37 8 568 9 99 10 851 11 382 ;
param: j :=  1 178 2 630 3 434 4 327 5 858 6 389 7 141 8 671 9 202 10 954 11 706 ;
param k0 :=  1 8  2 8  3 8  4 8  5 1  6 1  7 8  8 8  9 1  10 8 11 1;
param l0 :=  1 1  2 1  3 8  4 2  5 8  6 2  7 4  8 8  9 1  10 1 11 8;

# Execution, options
option solver ilogcp; solve;

# Post-Solution actions
display dist2;
display {i in IDX} (k[i],l[i]);

PS: If somebody asked me the meaning of a square dollar, I would not be able to answer!

JohnNichols · ‎12-18-2020

PS: If somebody asked me the meaning of a square dollar, I would not be able to answer!

Actually you would. If I define a meter it has length and zero thickness, if I define a square metre it has two lengths and the area if the two lengths are orthonormal is simple the length 1 times length 2, and the square indicates that we are using orthonormal measure, but money is just a counter, and if we square a counter clearly we get another count, but the count is the square and therefore is akin to the orthonormal with out the ortho or the normal. But we could define a $ and a UK pound, and multiple them

2 dollars and 2 pounds = 4 dollar pounds, and so we could in essence treat them as orthonormal.

Of course this works only because of the peculiar properties of the determinant.

And of course you were well aware of this problem, so I was wondering where did you get the idea.

but from somewhere

But just because we can do it mathematically in an abstract sense doesn't mean that it's useful in any way in the real world. Squared-dollars can't buy you lunch.

--------------------------------------------------------------------------------------

which is wrong because if I laid dollars as

A A

on a table, teh waiter would see 4 squares - made up of dollars and he/she would take them, placing them in a stack makes them 4 again.

The problem is mathematical equivalence.

mecej4 · ‎12-18-2020

JohnNichols:

Having said all that about square dollars, you seem to be quite ready for volume in space-time .

If you flip a $1 bill left->right, do you have a -$1 bill? If you rotate the $ bill 90 degrees, do have an imaginary dollar?

JohnNichols · ‎12-18-2020

!  Console1.f90 
!*********************************************************
!*                   SIMPLEX METHOD                      *
!*                   --------------                      *
!*                                                       *
!* LIST OF MAIN VARIABLES:                               *
!*                                                       *
!*  R:          MAXIMIZE = Y, MINIMIZE = N               *
!*  NV:         NUMBER OF VARIABLES OF ECONOMIC FUNCTION *
!*              (TO MAXIMIZE OR MINIMIZE).               *
!*  NC:         NUMBER OF CONSTRAINTS                    *
!*  TS:         SIMPLEX TABLE OF SIZE NC+1 x NV+1        *
!*  R1:         =1 TO MAXIMIZE, =-1 TO MINIMIZE          *
!*  R2:         AUXILIARY VARIABLE FOR INPUTS            *
!* NOPTIMAL     BOOLEAN IF FALSE, CONTINUE ITERATIONS    *
!*  XMAX:       STORES GREATER COEFFICIENT OF ECONOMIC   *
!*              FUNCTION.                                *
!*  RAP         STORES SMALLEST RATIO > 0                *
!*  V:          AUXILIARY VARIABLE                       *
!*  P1,P2:      LINE, COLUMN INDEX OF PIVOT              *
!*  XERR:       BOOLEAN IF TRUE, NO SOLUTION.            *
!* ----------------------------------------------------- *
!* PROBLEM DESCRIPTION:                                  *
!* A builder of houses can make 3 kinds of them with     *
!* various profits: 15000$, 17000$ and 20000$.           *
!* Each kind must respect following conditions:          *
!* 1) for supply reasons, the number of houses of kind 2 *
!*    built each month must not exceed the number of     *
!*    kind 3 by more than two.                           *
!* 2) for staff reasons, the buider can make each month  *
!*    up to 5, 5 and 3, respectively of kinds 1, 2 and 3.*
!* 3) for organisation reasons, the builder can make up  *
!*    to 8 houses monthly of kinds 1,2 and 3, respecti-  *
!*    vely in the proportions of 3, 2 and 1.             *
!* The builder, having these data, wants to maximize his *
!* monthly profit by knowing the number oh houses of     *
!* each kind to build monthly.                           *
!* ----------------------------------------------------- *
!* SAMPLE RUN:                                           *
!* (Maximize 15 X1 + 17 X2 + 20 X3 with conditions:      *
!*                    X2 -   X3 <= 2                     *
!*           3 X1 + 3 X2 + 5 X3 <= 15                    *
!*           3 X1 + 2 X2 +   X3 <= 8     )               *
!*                                                       *
!* LINEAR PROGRAMMING                                    *
!*                                                       *
!* MAXIMIZE ? Y                                          *
!*                                                       *
!* NUMBER OF VARIABLES OF ECONOMIC FUNCTION ? 3          *
!*                                                       *
!* NUMBER OF CONSTRAINTS ? 3                             *
!*                                                       *
!* INPUT COEFFICIENTS OF ECONOMIC FUNCTION:              *
!*       #1 ? 15                                         *
!*       #2 ? 17                                         *
!*       #3 ? 20                                         *
!*       Right hand side ? 0                             *
!*                                                       *
!* CONSTRAINT #1:                                        *
!*       #1 ? 0                                          *
!*       #2 ? 1                                          *
!*       #3 ? -1                                         *
!*       Right hand side ? 2                             *
!*                                                       *
!* CONSTRAINT #2:                                        *
!*       #1 ? 3                                          *
!*       #2 ? 3                                          *
!*       #3 ? 5                                          *
!*       Right hand side ? 15                            *
!*                                                       *
!* CONSTRAINT #3:                                        *
!*       #1 ? 3                                          *
!*       #2 ? 2                                          *
!*       #3 ? 1                                          *
!*       Right hand side ? 8                             *
!*                                                       *
!* RESULTS:                                              *
!*                                                       *
!*       VARIABLE #1:  0.333333                          *
!*       VARIABLE #2:  3.000000                          *
!*       VARIABLE #3:  1.000000                          *
!*                                                       *
!*       ECONOMIC FUNCTION: 76.000000                    *
!*                                                       *
!* (Building monthly 1/3, 3 and 1 house(s) of kinds 1,2, *
!*  3, the builder can make a monthly profit of 76000$). *
!* ----------------------------------------------------- *
!* REFERENCE:                                            *
!* Modèles pratiques de décision Tome 2, By Jean-Pierre  *
!* Blanger, PSI Editions, France, 1982.                  *
!*                                                       *
!*                F90 Release 1.0 By J-P Moreau, Paris.  *
!*                         (www.jpmoreau.fr)             *
!*********************************************************
PROGRAM SIMPLEX

implicit none
integer, parameter :: CMAX=10, VMAX=10

integer NC, NV, XERR;
real*8  TS(0:CMAX,0:VMAX)

  call Data(NV,NC,TS)
  call Simplex1(NV,NC,TS,XERR)
  call Results(NV,NC,TS,XERR)

  stop

END


  subroutine Data(NV,NC,TS)
    implicit none
    integer, parameter :: CMAX=10, VMAX=10
 !   real*8 TS(CMAX,VMAX)
    real*8  TS(0:CMAX,0:VMAX)
    real*8 R1,R2
    character*1 R
    integer NV, NC, J, I

	print *,' '
    print *,'LINEAR PROGRAMMING'
	print *,' '
    write(*,30,advance='no'); read *, R
	print *,' '
    write(*,40,advance='no'); read *, NV
	print *,' '
    write(*,50,advance='no'); read *, NC
    if (R.EQ.'Y'.OR.R.EQ.'y') then
      R1 = 1.0
    else
      R1 = -1.0
    end if
	print *,' '
    print *,'INPUT COEFFICIENTS OF ECONOMIC FUNCTION:'
    do J = 1, NV
	  write(*,10,advance='no') J; read *, R2
      TS(1,J+1) = R2 * R1
    end do
    write(*,60,advance='no'); read *, R2
    TS(1,1) = R2 * R1
    do I = 1, NC
      write(*,20) I
      do J = 1, NV
        write(*,10,advance='no') J; read *, R2
        TS(I + 1,J + 1) = -R2
      end do
      write(*,60,advance='no'); read *, TS(I+1,1)
    end do
	print *,' '
    print *,'RESULTS:'
	print *,' '
    do J=1, NV  
	  TS(0,J+1) = float(J)
    end do
    do I=NV+1, NV+NC
	  TS(I-NV+1,0) = float(I)
    end do
    return
10  format('       #',I1,' ? ')
20  format(' CONSTRAINT #',I1,' ? ')
30  format(' MAXIMIZE (Y/N) ? ')
40  format(' NUMBER OF VARIABLES OF ECONOMIC FUNCTION ? ')
50  format(' NUMBER OF CONSTRAINTS ? ')
60  format('       Right hand side ? ')
  end

  Subroutine Simplex1(NV,NC,TS,XERR)
  implicit none
    integer, parameter :: CMAX=10, VMAX=10
 !   real*8 TS(CMAX,VMAX)
    real*8  TS(0:CMAX,0:VMAX)
	integer P1,P2,XERR, NOPTIMAL
    
    integer NV, NC, J, I
    xerr = 0
    
10  call Pivot(NV,NC,TS,P1,P2)
    call Formula(NV,NC,TS,P1,P2)
    call Optimize(NV,NC,TS,NOPTIMAL,XERR)
    if (NOPTIMAL.EQ.1) goto 10
    return
  end

  Subroutine Pivot(NV,NC,TS,P1,P2)
  implicit none
    integer, parameter :: CMAX=10, VMAX=10
 !   real*8 TS(CMAX,VMAX)
    real*8  TS(0:CMAX,0:VMAX)
	integer P1,P2
    real*8 RAP,V,XMAX
    integer NV, NC, J, I
    
    XMAX = 0.d0
    do J=2, NV+1
	  if (TS(1,J) > 0.d0.AND.TS(1,J) > XMAX) then
        XMAX = TS(1,J)
        P2 = J
      end if
    end do
    RAP = 999999.d0
    do I=2, NC+1
      if (TS(I,P2) >= 0.d0) goto 10
      V = dabs(TS(I,1) / TS(I,P2))
      if (V < RAP) then
        RAP = V
        P1 = I
      end if
10  end do
    V = TS(0,P2)
	TS(0,P2) = TS(P1,0)
	TS(P1,0) = V
    return
  end

  Subroutine Formula(NV,NC,TS,P1,P2)
  implicit none
    integer, parameter :: CMAX=10, VMAX=10
 !   real*8 TS(CMAX,VMAX)
    real*8  TS(0:CMAX,0:VMAX)
    integer NV, NC, J, I
	integer P1,P2
    do I=1, NC+1
      if (I == P1) goto 70
      do J=1, NV+1
        if (J == P2) goto 60
        TS(I,J) = TS(I,J) - TS(P1,J) * TS(I,P2) / TS(P1,P2)
60    end do
70  end do
    TS(P1,P2) = 1.d0 / TS(P1,P2)
    do J=1, NV+1
      if (J == P2) goto 100
      TS(P1,J) = TS(P1,J) * dabs(TS(P1,P2))
100 end do
    do I=1, NC+1
      if (I == P1) goto 110
      TS(I,P2) = TS(I,P2) * TS(P1,P2)
110 end do
    return
  end   

  Subroutine Optimize(NV,NC,TS,NOPTIMAL,XERR)
  implicit none
    integer, parameter :: CMAX=10, VMAX=10
 !   real*8 TS(CMAX,VMAX)
    real*8  TS(0:CMAX,0:VMAX)
    integer NV, NC, J, I,NOPTIMAL
	integer XERR
    do I=2, NC+1
      if (TS(I,1) < 0.d0)  XERR = 1
    end do
    NOPTIMAL = 0
    if (XERR == 1) return
    do J=2, NV+1
      if (TS(1,J) > 0.d0)  NOPTIMAL = 1
    end do
    return
  end

  Subroutine Results(NV,NC,TS,XERR)
  implicit none
    integer, parameter :: CMAX=10, VMAX=10
    integer NV, NC, J, I
 !   real*8 TS(CMAX,VMAX)
    real*8  TS(0:CMAX,0:VMAX)
	integer XERR
    
!	if (XERR == 0) goto 30
!    print *,' NO SOLUTION.'!
!	goto 100
30  do I=1, NV
	  do J=2, NC+1
        if (TS(J,0).NE.float(I)) goto 70
        write(*,120) I, TS(J,1)
70    end do
    end do
    write(*,130) TS(1,1)
100 print *,' '
    return
120 format('       VARIABLE #',I1,': ',F10.6)
130 format('       ECONOMIC FUNCTION: ',F10.6) 
  end

! end of file simplex.f90

This is a simple simplex program, interestingly the original code would not run in any Fortran compiler I have ever seen. It took a while to get it fixed,

In terms of negative dollars,etc .. you are aware that -1 was a long argument in math, but we can define our new algebra.

I have no interest in space-time, I have a 13 year old daughter

JohnNichols · ‎12-18-2020

From the paper:

Let us consider a trivial physical system with Lb = 1.
Then the integral measures the volume of the corresponding piece of spacetime

As soon as math and physics people define a system and then simplify it. --- I love it.

JohnNichols · ‎12-18-2020

	integer XERR
    
!	if (XERR == 0) goto 30
!    print *,' NO SOLUTION.'!
!	goto 100
30    ...

I do not understand why if xerr came through as a zero this did not work. I kept getting the print no solution till I blanked it out.

mecej4 · ‎12-19-2020

That caution is borderline esoteric, and has to do with whether the boundaries of the feasible region are well-defined when inequality constraints are present. The question can be phrased "by how much does the l.h.s. have to be greater than the r.h.s. for the inequality to be considered as satisfied?"

You can make the caution go away by rephrasing your additional constraint as

subj to Limit {i in IDX}: 10*k[i]+l[i] >= m0[i]+1;

With this change, you should still see the same solution as before but without the caution.

Glad to see you rolling with AMPL!

JohnNichols · ‎12-19-2020

Many years ago I had to take a systems programming class. The professor was a really good guy, Dr Williams. He is the creator of the Filed Williams model, which is an early variant on Delft3d.

Anyway at the end of class he set the final assignment as four alternatives, 1 was a for a pass, 2 was for a credit 3 for a distinction and 4 for a high distinction. At my Uni the average GPA was 1.65 with a low standard deviation. My first year at Newcastle one young lady got straight HD's, she was the only one in 16000 students and everybody knew her name. She was a math major and had to be very bright.

No 4 was the travelling sales problem, so of course being an idiot I tried it. I could get it to solve with odd sized matrices but not evens. I played with it for a few weeks. By that time we had a computer for the students in the dept but we are shared a common folder for the class. One of the students, a good friend John Callaghan, wrote a simple program, we could only program in Fortran, had his working, so I opened his file, put a clear screen command and then a write statement saying "The Phantom Program snatcher strikes again." and then stop.

We are sitting around and John walks in, sits down and runs his program. I can still remember the swearing, I had to finally show him the problem. Of course at my current Uni I broke the honor code doing that and could be turfed out. Life was simpler then.

Dr Williams said to me, good solution, why did you not try the evens. But he still gave me a HD. So I was happy.

roy437 · ‎12-19-2020

Thank you mecej4,

Now I reformulated my problem and here is what the AMPL code looks like, which works very well.

# Adding transport in price, solving in integer numbers
# p[i] to be as close to pct[i] as possible

# Declarations
param M integer := 11;
param V integer := 579673;
set IDX = 1..M;
var p {IDX} integer >= 1;
param pft {IDX} integer >= 1;
param pct {IDX} integer >= 1;
param n {IDX} integer >= 1;

# Objective function
minimize dist2:
   sum {i in IDX} (p[i]-pct[i])^2;

# Constraints
subj to C: sum{i in IDX} n[i]*p[i] = V;
subj to Limit {i in IDX}: p[i] >= pft[i]+1;

#  Values of Parameters
data ;
param: n := 1 811 2 263 3 715 4 3 5 534 6 286 7 37 8 568 9 99 10 851 11 382;
param: pft := 1 22 2 240 3 61 4 10930 5 161 6 136 7 384 8 118 9 204 10 112 11 185;
param: pct := 1 23 2 253 3 64 4 11527 5 169 6 143 7 405 8 125 9 215 10 118 11 195;

# Execution, options
option solver ilogcp; solve;
display dist2;
display {i in IDX} p[i];

Solution:

ilogcp 12.10.0: ilogcp 12.10.0: optimal solution
49037 choice points, 45375 fails, objective 5
dist2 = 5

p[i] [*] :=
 1     23
 2    254
 3     63
 4  11527
 5    169
 6    143
 7    406
 8    125
 9    215
10    119
11    196
;

JohnNichols · ‎12-19-2020

What type of real problem are you solving?

roy437 · ‎12-20-2020

Here's the problem:
Calculate the new unit prices exactly to two decimal places so that the value of the shipping is distributed in proportion to the Val. In other words, the value of the shipping is included in the unit price of each product.

Solution: 55.12, 223.62, 603.25

andrew_4619 · ‎12-20-2020

I don't get that, the unit prices are a simple calculation that doesn't need all this complicated solving stuff or is this a different problem?

roy437 · ‎12-20-2020

Very often the result is with many decimals or even with an infinite number of decimals. It is desired to obtain those solutions that have only two decimals and that check the total invoice.

roy437 · ‎12-19-2020

Thanks John Nichols,

Here's another approach to the problem, "SIMPLEX METHOD".

roy437 · ‎12-19-2020

Thank you again mecej4

I introduced another parameter m0() and a constrain 10*k[i]+l[i] > m0[i] (did I write correctly?). I'm very happy that it works, but it gives me a message, what does that mean? Below you have the "mod" and the solution with "Caution:"

# Adding transport in price, solving in integer numbers
# {k, l} to be as close to {k0, l0} as possible

# Declarations
param M integer := 7;
param V integer := 398100;
set IDX = 1..M;
var k {IDX} integer >= 0 <= 9;
var l {IDX} integer >= 0 <= 9;
param k0 {IDX} integer;
param l0 {IDX} integer;
param n {IDX} integer >= 1;
param j {IDX} integer >= 0;
param m0 {IDX} integer >= 0;

# Objective function
minimize dist2:
   sum {i in IDX} ((k[i]-k0[i])^2 + (l[i]-l0[i])^2);

# Constraints
subj to C: sum{i in IDX} (100*j[i]+10*k[i]+l[i])*n[i] = V;
subj to Limit {i in IDX}: 10*k[i]+l[i] > m0[i];

#  Values of Parameters
data ;
param: n := 1 98 2 550 3 354 4 806 5 610 6 413 7 642;
param: j := 1 4 2 0 3 2 4 0 5 0 6 1 7 1;
param: k0 := 1 9 2 5 3 1 4 6 5 5 6 9 7 2;
param: l0 := 1 3 2 1 3 1 4 8 5 6 6 6 7 1;
param: m0 := 1 74 2 49 3 3 4 65 5 54 6 89 7 16;

# Execution, options
option solver ilogcp; solve;
display dist2;
display {i in IDX} (k[i],l[i]);

Solution:

price_with_transport.mod, line 22 (offset 607):
	Caution: Treating strict inequality constraint as a logical constraint.
context:  subj to Limit {i in IDX}: 10*k[i]+l[i] >  >>> m0[i]; <<< 

ilogcp 12.10.0: ilogcp 12.10.0: optimal solution
23558 choice points, 22247 fails, objective 6
dist2 = 6

: k[i] l[i]    :=
1   9    2
2   5    0
3   1    0
4   6    7
5   5    7
6   9    6
7   2    2
;

mecej4 · ‎12-17-2020

Arjen Markus is spot on, and I hope that his question gets you to think about methods other than blind brute force. The assignment appears to be a rather effective homework assignment, and it presents you with an opportunity to think/learn about the properties of modular arithmetic, and exploit that knowledge to make the program short, general and efficient.

You do not need the arbitrarily-deep DO loop capability for this problem. You should not use a CASE construct for m = 2, 3, .... Instead, make m a variable and write code that will work for any input m, perhaps up to a limit. Unless you have a valid reason to use floating point arithmetic, do not use it for a Diophantine problem.

Here is a hint to get you started. Find the answer by thinking and using pencil-and-paper:

What is the value of V mod 10 ?

If this is a homework assignment, you may be doing us as well as yourself a disservice by asking us to provide you with the code itself.

roy437 · ‎12-17-2020

This is not a homework assignment, it is a real problem that I face every day in my accounting work. Finding the solution with pencil and paper for m greater than 2 takes a long time.

JohnNichols · ‎12-18-2020

https://apps.dtic.mil/dtic/tr/fulltext/u2/a109487.pdf

Fortran Language