topic Re: Is it possible to redefine subroutines in a subroutine? in SAS/IML Software and Matrix Computations

Is it possible to redefine subroutines in a subroutine?

LauriN — Tue, 22 Jul 2014 11:23:17 GMT

I have this kind of external subroutine which uses original input vector p to make p-wl overlapping rolling windows of length wl. For each of those windows defined as "newp" I would like to use my subroutines. So they should occur after line " newp=p;" and p and m parameters in the subroutines should be newp and m1 respectfully as defined in this new subroutine "outer". Is it possible to for example make a module that compiles these subroutines and redefines their input parameters and then call this module in this "outer" subroutine and run this "mincvf" function? Logic of compiling module illustrated at the end:

start outer(p,e,pi,m); m=nrow(p);

wl=38; /window length/

m1=m-wl; /last window begins at m-wl/

newp=j(1,m1);

hyi=j(1,m1);

do x=1 to m1;

we=x+wl-1; /window end/

w=T(x:we);/window/

newp=p;

run comp (newp,m1,pi,e);

hyi=mincvf(CVF1);

end;

return (hyi);

finish;

The not-working compile module with the subroutines:

start comp(p,m,pi,e);
    start Kmod(x,h,pi,e);
    k=1/(h#(2#pi)##(1/2))#e##(-x##2/(2#h##2));
    return (k);
    finish;

start mhatx2(m,p,h,pi,e);

    t5=j(m,1);                   /mhatx omit x=t/
    do x=1 to nrow(p);
    i=T(1:m);
    temp1=x-i;
    ue=Kmod(temp1,h,pi,e)#p;
    le=Kmod(temp1,h,pi,e);
    t5=(sum(ue)-ue)/(sum(le)-le);
    end;
    return (t5);
    finish;

    Start CVF1(h) global(p,pi,e,m);
    cv3=j(m,1);
    cv3=1/nrow(p)#sum((p-mhatx2(m,p,h,pi,e))##2);
    return(cv3);
    finish;

    start mincvf(CVF1);
    optn={0,0};
    init=1;
    call nlpqn(rc, res,"CVF1",init);
    return (res);
    finish;
finish;

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Tue, 22 Jul 2014 14:30:30 GMT

In IML all subroutine definitions are global, even if you define then inside another subroutine. In your example, could you not accomplish the same goal by using a global clause in mhatx2 and Kmod? That is, use

start mhatx2(h) global(p,m,pi,e);

start Kmod(x,h) global(pi,e);

and manipulate p in your "outer" routine. maybe use something like:

start outer(p_input) global(p,e,pi,m); m=nrow(p_input);

[your existing code]

p = p_input(w);

hyi=mincvf(CVF1);

...

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Tue, 22 Jul 2014 14:36:34 GMT

Oops, I see that you will also have to manipulate m within the loop, but you get the idea?

Re: Is it possible to redefine subroutines in a subroutine?

LauriN — Wed, 23 Jul 2014 11:27:40 GMT

There were some mistakes like the subroutines were using m=nrow(p) instead of nrow(newp) as they now should and altogether the approach was somewhat complicated.

Here is what I did now, so trying to define the subroutines already outside to have those input values which "outer" function would define i.e. newp and nrow(newp), latter denoted now as m2. I also had such 2 lines to test the subroutines:

test1=T(1:38);

newp=p[test1];

and got the h for this first window as intended. When I ran the whole programm leaving those test value rows stay there as such, I actually got that m-wl times "NOTE: ABSGCONV convergence criterion satisfied." to my log. But ttt was m-wl times the constant h for the 1st window only. Strange.

I paste at the end the error message I get running the program as it's written below, doesn't include these 2 test lines.

proc iml;

EDIT kirjasto.pb var "p";

read all var "p" into p;

m=nrow(p);

m2=38;

pi=constant("pi");

e=constant("e");

start Kmod(x,h,pi,e);

k=1/(h#(2#pi)##(1/2))#e##(-x##2/(2#h##2));

return (k);

finish;

start mhatx2(m2,newp,h,pi,e);

t5=j(m2,1); /*mhatx omit x=t*/

do x=1 to m2;

i=T(1:m2);

temp1=x-i;

ue=Kmod(temp1,h,pi,e)#newp;

le=Kmod(temp1,h,pi,e);

t5=(sum(ue)-ue)/(sum(le)-le);

end;

return (t5);

finish;

Start CVF1(h) global(newp,pi,e,m2);

cv3=j(m2,1);

cv3=1/m2#sum((newp-mhatx2(m2,newp,h,pi,e))##2);

return(cv3);

finish;

start mincvf(CVF1);

optn={0,0};

init=1;

call nlpqn(rc, res,"CVF1",init);

return (res);

finish;

start outer(p,m) global(;

wl=38; /*window length*/

m1=m-wl; /*last window begins at m-wl*/

newp=j(wl,1);

hyi=j(m1,1);

do x=1 to m1;

we=x+wl-1; /*window end*/

w=T(x:we); /*window*/

newp=p;

hyi=mincvf(CVF1);

end;

return (hyi);

finish;

wl=38; /*window length*/

m1=m-wl; /*last window begins at m-wl*/

ttt=j(m1,1);

ttt=outer(p,m);

print ttt;

370 proc iml;

NOTE: IML Ready

371

372 EDIT kirjasto.basfraaka var "open";

373

374 read all var "open" into p;

375

376

377 m=nrow(p);

378 m2=38;

379 pi=constant("pi");

380 e=constant("e");

381

382

383 start Kmod(x,h,pi,e);

384 k=1/(h#(2#pi)##(1/2))#e##(-x##2/(2#h##2));

385 return (k);

386 finish;

NOTE: Module KMOD defined.

387

388

389 start mhatx2(m2,newp,h,pi,e);

390 t5=j(m2,1);

390! /*mhatx omit x=t*/

391 do x=1 to m2;

392 i=T(1:m2);

393 temp1=x-i;

394 ue=Kmod(temp1,h,pi,e)#newp;

395 le=Kmod(temp1,h,pi,e);

396 t5=(sum(ue)-ue)/(sum(le)-le);

397 end;

398 return (t5);

399 finish;

NOTE: Module MHATX2 defined.

400

401 Start CVF1(h) global(newp,pi,e,m2);

402 cv3=j(m2,1);

403 cv3=1/m2#sum((newp-mhatx2(m2,newp,h,pi,e))##2);

404 return(cv3);

405 finish;

NOTE: Module CVF1 defined.

406

407

408 start mincvf(CVF1);

409 optn={0,0};

410 init=1;

411 call nlpqn(rc, res,"CVF1",init);

412 return (res);

413 finish;

NOTE: Module MINCVF defined.

414

415

416 start outer(p,m);

417 wl=38;

417! /*window length*/

418 m1=m-wl;

418! /*last window begins at m-wl*/

419 newp=j(wl,1);

420 hyi=j(m1,1);

421 do x=1 to m1;

422 we=x+wl-1;

422! /*window end*/

423 w=T(x:we);

423! /*window*/

424 newp=p;

425 hyi=mincvf(CVF1);

426 end;

427 return (hyi);

428 finish;

NOTE: Module OUTER defined.

429

430

431 wl=38;

431! /*window length*/

432 m1=m-wl;

432! /*last window begins at m-wl*/

433

434 ttt=j(m1,1);

435 ttt=outer(p,m);

ERROR: (execution) Matrix has not been set to a value.

operation : [ at line 394 column 27

operands : newp, i

newp 0 row 0 col (type ?, size 0)

i 38 rows 1 col (numeric)

statement : START at line 383 column 1

traceback : module MHATX2 at line 383 column 1

module CVF1 at line 403 column 1

module MINCVF at line 408 column 1

module OUTER at line 425 column 1

ERROR: (execution) Matrix has not been set to a value.

operation : NLPQN at line 411 column 1

operands : *LIT1021, init

*LIT1021 1 row 1 col (character, size 4)

CVF1

init 1 row 1 col (numeric)

statement : CALL at line 411 column 1

traceback : module CVF1 at line 411 column 1

module MINCVF at line 408 column 1

module OUTER at line 425 column 1

436 print ttt;

ERROR: An exception has been encountered.

Please contact technical support and provide them with the following traceback information:

The SAS task name is [IML]

ERROR: Write Access Violation IML

Exception occurred at (01C89811)

Task Traceback

Address Frame (DBGHELP API Version 4.0 rev 5)

0000000001C89811 000000002403A880 tkmk:tkBoot+0x17AE1

0000000001C877BB 000000002403A8E0 tkmk:tkBoot+0x15A8B

0000000002F65527 000000002403A8E8 sashost:Main+0x1E097

0000000004014B05 000000002403A978 saswob:tkvercn1+0x3AC5

00000000079CB054 000000002403A9E8 saswobc:tkvercn1+0xA014

000000000401C748 000000002403AC90 saswob:tkvercn1+0xB708

0000000005BD4DCA 000000002403AC98 sasods:tkvercn1+0x1D3D8A

0000000005AF994D 000000002403B490 sasods:tkvercn1+0xF890D

0000000005BD4DCA 000000002403B4D0 sasods:tkvercn1+0x1D3D8A

0000000005A39BA7 000000002403B5C0 sasods:tkvercn1+0x38B67

0000000005BD4DCA 000000002403B600 sasods:tkvercn1+0x1D3D8A

0000000005AF9A82 000000002403BDC0 sasods:tkvercn1+0xF8A42

0000000005BD4DCA 000000002403BE00 sasods:tkvercn1+0x1D3D8A

0000000005AF36FF 000000002403C0B0 sasods:tkvercn1+0xF26BF

0000000005BD4DCA 000000002403C0F0 sasods:tkvercn1+0x1D3D8A

0000000005AF8378 000000002403C380 sasods:tkvercn1+0xF7338

0000000005BD4DCA 000000002403C3C0 sasods:tkvercn1+0x1D3D8A

0000000005AC0678 000000002403CAF0 sasods:tkvercn1+0xBF638

0000000005ABEDD0 000000002403CC80 sasods:tkvercn1+0xBDD90

0000000005BD4DCA 000000002403CCC0 sasods:tkvercn1+0x1D3D8A

0000000005A70257 000000002403D380 sasods:tkvercn1+0x6F217

0000000005BD4DCA 000000002403D3C0 sasods:tkvercn1+0x1D3D8A

0000000005AA7C04 000000002403DA30 sasods:tkvercn1+0xA6BC4

0000000023E938FA 000000002403DA38 sasiml:tkvercn1+0x628BA

0000000023E40DB0 000000002403E1B0 sasiml:tkvercn1+0xFD70

0000000023E3F8BD 000000002403F2A0 sasiml:tkvercn1+0xE87D

0000000023E3DE77 000000002403F4F0 sasiml:tkvercn1+0xCE37

0000000023E33587 000000002403FB90 sasiml:tkvercn1+0x2547

0000000023E312F1 000000002403FBF0 sasiml:tkvercn1+0x2B1

0000000002F5833B 000000002403FF20 sashost:Main+0x10EAB

0000000002F5DF7D 000000002403FF50 sashost:Main+0x16AED

00000000772859ED 000000002403FF58 kernel32:BaseThreadInitThunk+0xD

00000000773BC541 000000002403FF88 ntdll:RtlUserThreadStart+0x21

NOTE: The SAS System stopped processing this step because of errors.

NOTE: PROCEDURE IML used (Total process time):

real time 0.09 seconds

cpu time 0.07 seconds

Re: Is it possible to redefine subroutines in a subroutine?

Rick_SAS — Wed, 23 Jul 2014 12:31:27 GMT

I'd like to run your program, but you did not include the data. Can you either attach the kirjasto.basfraaka file or define the 'p' vector in the program?

Re: Is it possible to redefine subroutines in a subroutine?

LauriN — Wed, 23 Jul 2014 12:39:16 GMT

Sure, now there should be data attached!

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Wed, 23 Jul 2014 13:28:19 GMT

Your subroutine "outer" needs to declare newp in a global clause. Otherwise, when it sets newp it is working with a local variable only, and the "global" newp is not initialized. That is causing the first error in your log.

Re: Is it possible to redefine subroutines in a subroutine?

LauriN — Wed, 23 Jul 2014 13:45:58 GMT

Oh I see, how do I actually do that?

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Wed, 23 Jul 2014 14:50:55 GMT

start outer(p,m) global(newp);

Re: Is it possible to redefine subroutines in a subroutine?

LauriN — Wed, 23 Jul 2014 15:17:06 GMT

Oh it works "that way" also, good to know! Thanks. Now it's producing at least something. The program also gives now such a warning, few log-fulls of them. Seems to refer to mhatx2 subroutine and to (sum(le)-le) but it shouldn't actually be zero at any point. Kmod doesn't go to zero at any input x-i.

start mhatx2(m2,newp,h,pi,e);

t5=j(m2,1); /*mhatx omit x=t*/

do x=1 to m2;

i=T(1:m2);

temp1=x-i;

ue=Kmod(temp1,h,pi,e)#newp;

le=Kmod(temp1,h,pi,e);

t5=(sum(ue)-ue)/(sum(le)-le);

end;

return (t5);

finish;

WARNING: Division by zero, result set to missing value.

operation : / at line 463 column 22

operands : _TEM1003, _TEM1006

_TEM1003 1 row 1 col (numeric)

_TEM1006 1 row 1 col (numeric)

statement : ASSIGN at line 463 column 1

traceback : module MHATX2 at line 463 column 1

module CVF1 at line 470 column 1

module MINCVF at line 475 column 1

module OUTER at line 492 column 1

WARNING: Division by zero, result set to missing value.

count : number of occurrences is 2

operation : / at line 463 column 22

operands : _TEM1003, _TEM1006

_TEM1003 1 row 1 col (numeric)

_TEM1006 1 row 1 col (numeric)

statement : ASSIGN at line 463 column 1

traceback : module MHATX2 at line 463 column 1

module CVF1 at line 470 column 1

module MINCVF at line 475 column 1

module OUTER at line 492 column 1

WARNING: Division by zero, result set to missing value.

operation : / at line 463 column 22

operands : _TEM1003, _TEM1006

_TEM1003 1 row 1 col (numeric)

_TEM1006 1 row 1 col (numeric)

statement : ASSIGN at line 463 column 1

traceback : module MHATX2 at line 463 column 1

module CVF1 at line 470 column 1

module MINCVF at line 475 column 1

module OUTER at line 492 column 1

WARNING: Division by zero, result set to missing value.

count : number of occurrences is 2

operation : / at line 463 column 22

operands : _TEM1003, _TEM1006

_TEM1003 1 row 1 col (numeric)

_TEM1006 1 row 1 col (numeric)

statement : ASSIGN at line 463 column 1

traceback : module MHATX2 at line 463 column 1

module CVF1 at line 470 column 1

module MINCVF at line 475 column 1

module OUTER at line 492 column 1

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Wed, 23 Jul 2014 19:55:16 GMT

the warning appears to be saying that the value of (sum(le)-le) is 0. This could be due to loss

of precision if le >> le, i != x. I recommend you add some conditional print statements, like

if (sum(le) - le) = 0 then do;

print temp1, x, i, ue, le;

end;

in front of the statement: t5=(sum(ue)-ue)/(sum(le)-le);

to see what is going on.

Re: Is it possible to redefine subroutines in a subroutine?

LauriN — Thu, 24 Jul 2014 11:54:21 GMT

Hmm interesting. x takes normally values between 1 and 38 and i is also 1 to 38 and thereby temp1 takes values between -37 and 38.. Anyway then the problem apparently is in minimizing function because if h approximates 0, then Kmod does so also. Have to apply restriction to it. How might I limit the desirable solutions?

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Thu, 24 Jul 2014 20:22:38 GMT

You could specify a lower bound or linear constraint, see in the IML doc under Nonlinear Optimization Examples->Details->Parameter Constraints

Re: Is it possible to redefine subroutines in a subroutine?

LauriN — Fri, 25 Jul 2014 09:51:42 GMT

So because h=0.4, "le" is some 0.08 and decreases steeply going lower in h, I use:

conh={0.4 . .,. . .,. . .};

start mincvf(CVF1);

optn={0,0};

init=1;

call nlpqn(rc, res,"CVF1",init) blc="conh";

return (res);

finish;

Is it missing something? The divided by zero error keeps happening and "outer" produces values under 0.4 still which it shouldn't.

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Fri, 25 Jul 2014 15:43:54 GMT

The optimization method should only be generating feasible points in the search. It looks like you are setting conh before you define mincvf(); your statement:

conh={0.4 . .,. . .,. . .};

should come AFTER the"start mincvf(CVF1);" statement.

Re: Is it possible to redefine subroutines in a subroutine?

LauriN — Sun, 27 Jul 2014 12:40:57 GMT

This way it works, had to put that option and constrain vector both into the input parentheses:

Now I get no division by 0 warning. The previously miss-specified-due-loss-of-precision point's are now not specified at all and the value is substituted by 0.14 but the error isn't likely big.

start mincvf(CVF1);

con={0.14 . .,. . .,. . .};

optn={0,0};

init=1;

call nlpqn(rc, res,"CVF1",init,optn,con);

return (res);

finish;

Thanks for the help!