SandMatrix Released - The trilogy is complete.



#2

The trilogy of Math modules is finally complete - including manuals and even real Overlays... (some still available)

- 41Z

- SandMath

- SandMatrix Manual

Enjoy!

PS.- Glad to see the 41 platform in such a good health.... now it's time for an extended break!

Edited: 28 Aug 2013, 12:29 p.m.


#3

Hey Angel. Thanks again for all you do.

I'd like to make a request / suggestion. :-)

How about showing us what a linear programming solution would look like using the SandMatrix module?

Solving something like:

Max 2x+4y
Subject to:
X + 3Y <= 5000
Y > 100
5X + Y <= 3000
X > 0
Y > 0

or such. You know, a classical optimization situation. Mix or max. I would think the coefficients could be stored in a matrix and a program given the mix or max choice and the number of variables, etc.

Why you?

Because you are da man. lol. :-)


#4

Sorry to disappoint you Gene but - I really think this assignment should be given to a more capable person. I wouldn't touch the linear optimization with a 12-foot pole, sort of always hated the SIMPLEX method and its derivatives :(

Suffice it to say some things we're good at and some others, err... you get the picture. Besides I'm in my break now, right?

Best,
'AM

#5

Gene,

A little off-topic, but I'm trying to avoid real work.

A while back I wrote an implementation of the downhill simplex on the WP-34s. While that is generally used for non-linear optimization and has a host of caveats so it takes some massaging of the problem to get it to work, I nonetheless thought that I would give it a try with the example in your posting. It did give (within the tolerance value) the same solution as WolframAlpha (http://goo.gl/DP2BU). It took about 1000 iterations of the simplex. I ran it on the emulator so I am not sure how fast it is on the actual hardware.

I know that your question was about SandMatrix, but I thought that as a curiosity I would toss out my solution. To run it, you XEQ'NMS' with the number of variables (2) in X on the stack.

To use it, you have to code up a function called FUN which implements your function along with constraints. Because my implementation of the simplex assumes minimization, I changed the function value to be negated. To handle constraints in your example, I simply added a large "penalty" to the value of the function if a constraint is not met.

The NMS routine is not really quite ready to be called general purpose. It currently does not let you enter a starting point (it starts around 0,0) and it does not let you specify the tolerance which is currently hardwired. As soon as I clean those up one of these days I'll probably post it to the articles section.

However, the cool thing is that you can use it to minimize pretty much any type of function(s) with up to 8 variables. (All the traditional pitfalls of Nelder-Mead apply).

The listings below are meant to be run through the wp34s pre-processor / assembler.

Here is the FUN function that codes up your example:

LBL'FUN'
LocR 4

STO .00
R[v]
STO .01

/* big value for constraint hack */
9
EEX
9
STO .03

/* evaluate 2x+4y but negate it so that we can minimize */
2
RCL* .00
4
RCL* .01
+
CHS
STO .02

/* x+3y<=5000 */
5
0
0
0
ENTER
3
RCL* .01
RCL+ .00
x>? Y
XEQ 98

/* y>100 */
# 100
RCL .01
x<=? Y
XEQ 98

/* 5x+y<=3000 */
3
0
0
0
ENTER
5
RCL* .00
RCL+ .01
x>? Y
XEQ 98

/* x > 0 */
0
RCL .00
x<=? Y
XEQ 98

/* y > 0 */
0
RCL .01
x<=? Y
XEQ 98

RCL .02
RTN

/* add error for unmet constraint*/
LBL 98
RCL .03
STO+ .02
RTN

And, here is the routine that does the solution. The comments make sense if you follow along with the amoeba code in Numerical Recipes.

/* 
Nelder Mead (Downhill Simplex) Method
Multi-dimensional Minimization

Limitations: The code supports up to 8 dimensions, but you will
probably not be able to use that many because of
memory restrictions. It depends on the size of your
evaluation function, register allocations, etc.

Global Register required: N^2 + N + 20

Stack: Requires a stack size of 8

Local Registers: (19 at deepest point) + (your eval function)

Prerequisite: You need to have a global function called 'FUN'
which evaluates your function. The parameters to FUN
are all passed on the stack. FUN needs to return
one value (the evaluated function) in X.

Input: (X) Number of dimensions

Output: R03: F(x1, x2, x3....) at minimum
R04: x1
R05: x2
R06: X3
...

Configuration: This algorithm is known for getting stuck, especially
if the function has lots of local minima or
discontinuities. There is a configuration value that
defines the maximum number of times it shall adjust
the simplex before giving up. It defaults to 10000.
If you wish to adjust this value, it is assigned to
local register .06 right at the beginning of NMS.

The initial guess is aribrarily picked to be between
0 and 1 on all dimensions. If you wish to change the
initial values, change the initValues function
at LBL 31. Make sure that when that function returns
it has valid values in the Y array for the
simplex vertices you have chosen as the starting point.

Algorithm: This is a direct derivative of the Downhill Simplex
method as described and implemented in Numerical
Recipes. Please see that text for a description
of how the algorithm works.
*/

LBL'NMS'
LocR 007
REGS 100
SSIZE8
CLREGS

/* GREG[GREG_DIM] = NP; */
/* the number of dimensions got passed in in the X register */
STO 00

/* init NMAX */
1
SDL 4
STO .06

/* initValues() */
XEQ 31

/* calcPsum(); */
XEQ 30

/* loop indefintely until exit condition */
/* for (;;) { */
mainLoop:: 0

/* LREG[LREG_AMOEBA_ILO] = 0; */
STO .02

/* GREG[GREG_IHI] = getY(0) >
getY(1) ? (LREG[LREG_AMOEBA_INHI]=1,0) (LREG[LREG_AMOEBA_INHI]=0,1); */
0
XEQ 20
1
XEQ 20
x<? Y
JMP m1
0
STO .03
1
STO 02
JMP m2
m1:: 1
STO .03
0
STO 02

/* for (LREG[LREG_AMOEBA_I]=0;
LREG[LREG_AMOEBA_I] < (int)GREG[GREG_DIM]+1;
LREG[LREG_AMOEBA_I]++) { */
m2:: RCL 00
STO .00

mainLoop1:: RCL .00

/* if (getY((int)LREG[LREG_AMOEBA_I]) <= getY((int)LREG[LREG_AMOEBA_ILO])) */
XEQ 20
RCL .02
XEQ 20
x<? Y
JMP m3

/* LREG[LREG_AMOEBA_ILO] = (int)LREG[LREG_AMOEBA_I]; */
RCL .00
STO .02

/* if (getY((int)LREG[LREG_AMOEBA_I]) > getY((int)GREG[GREG_IHI])) { */
m3:: RCL .00
XEQ 20
RCL 02
XEQ 20
x>=? Y
JMP m4

/* LREG[LREG_AMOEBA_INHI]=(int)GREG[GREG_IHI]; */
RCL 02
STO .03

/* GREG[GREG_IHI]=LREG[LREG_AMOEBA_I]; */
RCL .00
STO 02
JMP m5

/* else if (getY((int)LREG[LREG_AMOEBA_I]) >
getY((int)LREG[LREG_AMOEBA_INHI]) &&
LREG[LREG_AMOEBA_I] != GREG[GREG_IHI]) & */
m4:: RCL .00
XEQ 20
RCL .03
XEQ 20
x>=? Y
JMP m5
RCL 02
RCL .00
x=? Y
JMP m5

/* LREG[LREG_AMOEBA_INHI] = (int)LREG[LREG_AMOEBA_I]; */
STO .03

/* loop */
m5:: DSL .00
JMP mainLoop1

/* If tolerance is met, put the appropriate values into the output and be gone */
/* if (FTOL > 2.0*fabs(getY((int)GREG[02])-
getY((int)LREG[02]))/
(fabs(getY((int)GREG[02]))+fabs(getY((int)LREG[02]))+TINY)) { */

RCL 02
XEQ 20
RCL .02
XEQ 20
-
ABS
RCL 02
XEQ 20
ABS
RCL .02
XEQ 20
ABS
+
1
SDR 10
+
/
2
*

1
EEX
CHS
1
6
x<=? Y
JMP m6

/* GREG[GREG_SOLUTION_Y] = getY(LREG[LREG_AMOEBA_ILO]) */
RCL .02
XEQ 20
STO 03

/* for (LREG[LREG_AMOEBA_I]=0; LREG[LREG_AMOEBA_I] < GREG[GREG_DIM]; LREG[LREG_AMOEBA_I]++) { */
RCL 00
STO .00
DEC .00

/* GREG[04 + (int)LREG[00]] = getP((int)LREG[02], (int)LREG[00]); */
mainLoop2:: RCL .02
RCL .00
XEQ 24
4
RCL+ .00
x[<->] Y
STO[->]Y

/* loop */
DSL .00
JMP mainLoop2

/* done */
GTO 99

/* if ((int)GREG[GREG_NUMCALLS] >= NMAX) */
m6:: RCL 01
x<? .06
JMP m7
CL[alpha]
"Loop Max"
VIEW[alpha]
GTO 99

/* GREG[GREG_NUMCALLS] += 2; */
m7:: INC 01
INC 01

/* LREG[LREG_AMOEBA_YTRY]=amotry(-1.0); */
1
CHS
XEQ 26
STO .05

/* if (LREG[LREG_AMOEBA_YTRY] <= getY((int)LREG[LREG_AMOEBA_ILO])) */
RCL .02
XEQ 20
x<? Y
JMP m8
2
XEQ 26
STO .05
JMP m9

/*if (LREG[LREG_AMOEBA_YTRY] >= getY((int)LREG[LREG_AMOEBA_INHI])) { */
m8:: RCL .03
XEQ 20
RCL .05
x<? Y
JMP m10

/* LREG[LREG_AMOEBA_YSAVE] = getY((int)GREG[GREG_IHI]); */
RCL 02
XEQ 20
STO .04

/* LREG[LREG_AMOEBA_YTRY]=amotry(0.5); */
# 1/2
XEQ 26
STO .05

/* if (LREG[LREG_AMOEBA_YTRY] >= LREG[LREG_AMOEBA_YSAVE]) { */
RCL .04
x>? Y
JMP m9

/* for (LREG[LREG_AMOEBA_I]=0; LREG[LREG_AMOEBA_I] < (int)GREG[GREG_DIM]+1 ;LREG[LREG_AMOEBA_I]++) { */
RCL 00
STO .00
mainLoop3:: RCL .00

/* if (LREG[LREG_AMOEBA_I] != (int)LREG[LREG_AMOEBA_ILO]) { */
x=? .02
JMP m11

/* for (LREG[LREG_AMOEBA_J] = 0; LREG[LREG_AMOEBA_J] < GREG[GREG_DIM]; LREG[LREG_AMOEBA_J]++) { */
RCL 00
STO .01
DEC .01
mainLoop4:: RCL .00

/*setPsum((int)LREG[LREG_AMOEBA_J], 0.5 *
(getP((int)LREG[LREG_AMOEBA_I], (int)LREG[LREG_AMOEBA_J]) +
getP((int)LREG[LREG_AMOEBA_ILO], (int)LREG[LREG_AMOEBA_J]))); */
RCL .01
XEQ 24
RCL .02
RCL .01
XEQ 24
+
2
/
RCL .01
x[<->] Y
XEQ 23

/* setP((int)LREG[LREG_AMOEBA_I], (int)LREG[LREG_AMOEBA_J], getPsum((int)LREG[LREG_AMOEBA_J])); */
RCL .00
RCL .01
RCL .01
XEQ 22
XEQ 25

/* increment and loop */
DSL .01
JMP mainLoop4

/* setY((int)LREG[LREG_AMOEBA_I], func(&GREG[(int)GREG_PSUM])); */
RCLS 12
XEQ'FUN'
RCL .00
x[<->] Y
XEQ 21

/* increment and loop */
m11:: DSL .00
JMP mainLoop3

/* GREG[GREG_NUMCALLS] += GREG[GREG_DIM]; */
RCL 00
STO+ 01
XEQ 30

JMP m9

m10:: DEC 01

/* back to the top of the loop */
m9:: JMP mainLoop

/* we are done */
LBL 99
STOP
RTN

/* get the value pointed to by X, and lose the pointer from the stack */
LBL 00
RCL[->]X
x[<->] Y
DROP
RTN

/* getY */
LBL 20
3
+
XEQ 00
RTN

/* setY */
LBL 21
x[<->] Y
3
+
x[<->] Y
STO[->]Y

RTN

/* getPsum */
LBL 22
# 12
+
XEQ 00
RTN

/* setPsum */
LBL 23
x[<->] Y
# 12
+
x[<->] Y
STO[->]Y
RTN

/* getP */
LBL 24
x[<->] Y
RCL[times] 00
+
# 20
+
XEQ 00
RTN

/* setP */
LBL 25
[<->] ZYXT
RCL[times] 00
+
# 20
+
x[<->] Y
STO[->]Y
RTN

/* amotry */
LBL 26
LocR 13 /* max NP of 8 */

/* LREG[LREG_AMOTRY_IN_FAC] = fac */
STO .04

/* LREG[LREG_AMOTRY_FAC1] = (1.0 - LREG[LREG_AMOTRY_IN_FAC]) / GREG[GREG_DIM]; */
CHS
1
+
RCL 00
/
STO .01

/* LREG[LREG_AMOTRY_FAC2] = LREG[LREG_AMOTRY_FAC1]-LREG[LREG_AMOTRY_IN_FAC]; */
RCL .04
-
STO .02

/* for (LREG[LREG_AMOTRY_J] = 0; LREG[LREG_AMOTRY_J] < GREG[GREG_DIM]; LREG[LREG_AMOTRY_J]++) { {*/
RCL 00
STO .00
DEC .00

/* LREG[(int)(LREG_AMOTRY_PTRYOFFSET +
LREG[LREG_AMOTRY_J])] = getPsum((int)LREG[LREG_AMOTRY_J]) *
LREG[LREG_AMOTRY_FAC1] -
getP((int)GREG[GREG_IHI], (int)LREG[LREG_AMOTRY_J]) * LREG[LREG_AMOTRY_FAC2]; */
amotryLoop1:: RCL .00
XEQ 22
RCL .01
[times]
RCL 02
RCL .00
XEQ 24
RCL .02
[times]
-
RCL .00
# 117
+
x[<->] Y
STO[->]Y

/* increment the loop variable and loop */
DSL .00
JMP amotryLoop1

/*LREG[LREG_AMOTRY_YTRY] = func(&LREG[LREG_AMOTRY_PTRYOFFSET]); */
RCLS .05
XEQ'FUN'
STO .03

/* if (LREG[LREG_AMOTRY_YTRY] < getY((int)GREG[GREG_IHI])) { { */
RCL 02
XEQ 20
RCL .03
x>=? Y
JMP amotrycond1

/* setY((int)GREG[GREG_IHI], LREG[LREG_AMOTRY_YTRY]); */
RCL 02
RCL .03
XEQ 21

/* for (LREG[LREG_AMOTRY_J] = 0; LREG[LREG_AMOTRY_J] < GREG[GREG_DIM]; LREG[LREG_AMOTRY_J]++) { { */
RCL 00
STO .00
DEC .00

/* setPsum((int)LREG[LREG_AMOTRY_J], getPsum((int)LREG[LREG_AMOTRY_J]) + LREG[(int)
(LREG_AMOTRY_PTRYOFFSET + LREG[LREG_AMOTRY_J])] - getP((int)GREG[GREG_IHI],
(int)LREG[LREG_AMOTRY_J])); */
amotryLoop2:: RCL .00
ENTER

/* setPsum((int)LREG[LREG_AMOTRY_J], getPsum((int)LREG[LREG_AMOTRY_J]) +
LREG[(int)(LREG_AMOTRY_PTRYOFFSET + LREG[LREG_AMOTRY_J])] -
getP((int)GREG[GREG_IHI], (int)LREG[LREG_AMOTRY_J])); */
XEQ 22
# 117
RCL+ .00
XEQ 00
+

RCL 02
RCL .00
XEQ 24
-

XEQ 23

/* setP((int)GREG[GREG_IHI], (int)LREG[LREG_AMOTRY_J], LREG[(int)(LREG_AMOTRY_PTRYOFFSET + LREG[LREG_AMOTRY_J])]); */
RCL 02
RCL .00

# 117
RCL+ .00
XEQ 00
XEQ 25

/* increment the loop variable and loop */
DSL .00
JMP amotryLoop2

/* return LREG[LREG_AMOTRY_YTRY]; */
amotrycond1:: RCL .03
RTN

/* calcPsum */
LBL 30
LocR 3

/* for (LREG[LREG_CALCSUM_J] = 0; LREG[LREG_CALCSUM_J] < GREG[GREG_DIM]; LREG[LREG_CALCSUM_J]++) { */
RCL 00
STO .01
DEC .01
calcSumLoop1:: 0
/* for (LREG[LREG_CALCSUM_SUM] = 0.0, LREG[LREG_CALCSUM_I]=0; LREG[LREG_CALCSUM_I] <
(int)GREG[GREG_DIM]+1; LREG[LREG_CALCSUM_I]++) */
STO .02
RCL 00
STO .00
calcSumLoop2:: RCL .00

/* LREG[LREG_CALCSUM_SUM] = LREG[LREG_CALCSUM_SUM] + getP((int)LREG[LREG_CALCSUM_I], (int)LREG[LREG_CALCSUM_J]); */
RCL .01
XEQ 24
STO+ .02

/* increment loop 2 and loop */
DSL .00
JMP calcSumLoop2

/* setPsum((int)LREG[LREG_CALCSUM_J], LREG[LREG_CALCSUM_SUM]); */
RCL .01
RCL .02
XEQ 23

/* loop */
DSL .01
JMP calcSumLoop1

RTN

/* initValues */
LBL 31
LocR 10

/* for (LREG[LREG_INIT_I]=0; LREG[LREG_INIT_I] < (GREG[GREG_DIM]+1); LREG[LREG_INIT_I]++) { */
RCL 00
STO .00
initValuesLoop1:: RCL 00
STO .01
DEC .01
/* for (LREG[LREG_INIT_J] = 0;LREG[LREG_INIT_J] < GREG[GREG_DIM]; LREG[LREG_INIT_J]++) { { */
initValuesLoop2:: RCL .01
/* LREG[LREG_INIT_XOFFSET + (int)LREG[LREG_INIT_J]] = (LREG[LREG_INIT_I] == (LREG[LREG_INIT_J]+1) ? 1.0 : 0.0); */
0
+
RCL .00
0
x[<->] Y
x=? Z
INC Y
DROP

RCL .01
# 114
+
x[<->] Y
STO[->]Y

/* setP((int)LREG[LREG_INIT_I], (int)LREG[LREG_INIT_J], LREG[LREG_INIT_XOFFSET + (int)LREG[LREG_INIT_J]]); */
RCL .00
RCL .01
[<->] ZXYT
XEQ 25

/* increment loop2 and loop */
DSL .01
JMP initValuesLoop2

/* setY((int)LREG[LREG_INIT_I],func(&LREG[LREG_INIT_XOFFSET])); */
RCLS .02

XEQ'FUN'
RCL .00
x[<->] Y
XEQ 21

/* increment loop1 and loop */
DSL .00
JMP initValuesLoop1

RTN

Edited: 30 Aug 2013, 5:46 p.m.


Possibly Related Threads...
Thread Author Replies Views Last Post
  Droid48 Reader 1.2 Released Chris Dreher 0 141 10-23-2013, 02:26 PM
Last Post: Chris Dreher
  Latest & Last SandMath is released Ángel Martin 0 175 07-02-2013, 12:44 AM
Last Post: Ángel Martin
  The 'Complete' Collection Keith Midson 8 699 05-12-2013, 09:52 PM
Last Post: Namir
  HP 50g: Major updates for MLP / OSE / HLP released Software49g 2 203 04-15-2013, 04:49 AM
Last Post: Michael Lopez
  My Mostek collection is complete Michael de Estrada 15 963 04-10-2013, 10:01 PM
Last Post: Michael de Estrada
  Official released of the TI-84 Plus C Silver Edition Mic 12 523 01-21-2013, 01:38 PM
Last Post: Eric Smith
  Droid48 Reader 1.1 Released Chris Dreher 0 131 01-15-2013, 12:43 AM
Last Post: Chris Dreher
  Complete list of PPC members Gonzalo Fernandez 0 141 01-04-2013, 03:01 PM
Last Post: Gonzalo Fernandez
  SpeedUI: Version 12.05 Holiday Edition released:-) Raymond Del Tondo 3 236 12-30-2012, 02:18 AM
Last Post: Han
  WP 34S Owner's Manual released for print Walter B 28 1,030 12-10-2012, 08:36 AM
Last Post: Walter B

Forum Jump: