Is it possible to redefine subroutines in a subroutine?

[LauriN]

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;

[Hutch_sas]

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);

...

[Hutch_sas]

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

[LauriN]

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)

         1

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

[Rick_SAS]

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?

[LauriN]

Sure, now there should be data attached!

[Hutch_sas]

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.

[LauriN]

Oh I see, how do I actually do that?

[Hutch_sas]

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

[LauriN]

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)

         0

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)

         0

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)

         0

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)

         0

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

[Hutch_sas]

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.

[LauriN]

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?

[Hutch_sas]

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

[LauriN]

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.

[Hutch_sas]

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.