home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / sa104os2.zip / SATHR104.ZIP / SATHER / LIBRARY / MAT.SA < prev    next >
Text File  |  1995-02-04  |  16KB  |  457 lines
  1. -- Copyright (C) International Computer Science Institute, 1994.  COPYRIGHT  --
  2. -- NOTICE: This code is provided "AS IS" WITHOUT ANY WARRANTY and is subject --
  3. -- to the terms of the SATHER LIBRARY GENERAL PUBLIC LICENSE contained in    --
  4. -- the file "Doc/License" of the Sather distribution.  The license is also   --
  5. -- available from ICSI, 1947 Center St., Suite 600, Berkeley CA 94704, USA.  --
  6. --------> Please email comments to "sather-bugs@icsi.berkeley.edu". <----------
  7.  
  8. -- mat.sa: Matrices of FLT.
  9. -------------------------------------------------------------------
  10. class MAT is
  11.    -- Matrices of FLT.
  12.    
  13.    include ARRAY2{FLT};
  14.       -- Gets 
  15.       -- create(nr,nc),create(ARRAY{ARRAY{T}}), 
  16.       -- copy
  17.       -- aget, aset, 
  18.       -- nr, nc
  19.       -- row_ind!, col_ind!
  20.       -- row_elt!, col_elt!
  21.       -- set_row!, set_col!
  22.       -- diag_elt! set_diag_elt!
  23.       -- to_portion_of
  24.    
  25.    create_identity(sz:INT):SAME is
  26.       -- A new square identity matrix of size sz.
  27.       return(create_diagonal(sz,1.0)) end;
  28.  
  29.    create_diagonal(sz: INT, diag_val: FLT): SAME
  30.       -- A diagonal matrix with diag_val on the diagonal
  31.       pre sz>0 is
  32.       res ::= #SAME(sz,sz); 
  33.       loop i::=sz.times!; res[i,i]:=diag_val end; 
  34.       return(res);
  35.       end;
  36.    
  37.    same_shape_p(m1:SAME):BOOL is
  38.       -- True if self and m1 have the same dimensions.
  39.       return(m1.nr=nr and m1.nc=nc) end;
  40.    
  41. --   row_col_ind!: TUP{INT,INT} is
  42. --      loop r::=asize1.dotimes!;
  43. --     loop yield(#TUP2{INT,INT}(r,asize2.dotimes!)) end; end; end;
  44.    
  45.    to(val:SAME) 
  46.       -- Set self equal to val
  47.       pre same_shape_p(val) is
  48.       loop set!(val.elt!) end end;
  49.  
  50.    row_elts_in(r:INT, s:FSTR): FSTR is
  51.       -- Write successive elements of row `r' into s.
  52.       res ::= s;
  53. --    loop res := res + row_elt!(r).str end end;                                -- NLP ?
  54.       loop res := res + row_elt!(r).str end; return res; end;                   -- NLP ?
  55.  
  56. --   row_lists_in(s:FSTR) is
  57. --      -- Write lists containing successive rows into s.
  58. --      loop s.listify(ITER(row_elts_in(nr.do,_))) end end;
  59.  
  60. --   str_in(STRCAT) is
  61. --      -- Write a list version of self into arg.
  62. --      arg.listify(ITER(row_lists_in(_))) end;
  63.  
  64. --   str:STR is
  65. --      -- A string version of self.
  66. --      return(#FSTR+str; (sc); res:=sc.str end;
  67.  
  68. --   append_row_lists!(STRCAT) is
  69. --      -- Appends list form of successive rows to arg.
  70. --      loop
  71. --      end end;
  72.  
  73.    
  74.    str:STR is
  75.       -- The string form of self represented as a list of rows,
  76.       -- eg. "||1.00,2.33|,|4.5,2.8|,|9.7,3.2||".
  77.       sc::=#FSTR + "|";
  78.       loop r::=nr.times!; 
  79.      if r/=0 then sc := sc+ "," end; 
  80.      sc := sc + "|";
  81.      loop c::=nc.times!; 
  82.         if c/=0 then sc := sc+"," end; 
  83.         sc := sc + [r,c]; end;
  84. --        #OUT+r+","+c+":"+sc+"\n"; end;
  85.      sc := sc+ "|" end; 
  86.       sc := sc+"|"; 
  87.       return(sc.str) end;
  88.       
  89. --   from_str(STR):SAME is
  90.       -- Create a new matrix from a string with a list of rows,
  91.       -- eg. "((1.00 2.33) (4.5 2.8) (9.7 3.2))".
  92.       -- Returns `void' if not of this form (extra white space is ignored).
  93.       -- Must be at least one entry.
  94. --      sc::=arg.cursor; tmp::=LIST{LIST{FLT}}::create;
  95. --      sc.skip_space; if not sc.c_accept('(') then return end;
  96. --      loop sc.skip_space; if sc.c_accept(')') then break end;
  97. --     tmp:=tmp.push(LIST{FLT}::create);
  98. --     if not sc.c_accept('(') then return end;
  99. --     loop sc.skip_space; if sc.c_accept(')') then break end;
  100. --        tmp:=tmp.push(tmp.pop.push(sc.get_f));
  101. --        if sc.error/=0 then return end end end;
  102. --      loop if tmp.elts.size/=tmp[0].size then return end;
  103. --      res:=create(tmp.size,tmp.top.size);
  104. --      loop r::=res.nr.do; loop res.set_row_elts(r,tmp[r].elts) end end end;
  105.         
  106. --   to_portion_of(SAME) is
  107. --      -- Copy into self as much of arg as will fit and return it. Don't
  108. --      -- alter other elements.
  109. --      loop r::=nr.min(arg.nr).dotimes;
  110. --     loop set_row_elts(r,arg.row_elts) end end end;
  111.    
  112.    scale_by(s:FLT) is
  113.       -- Scale self by s.
  114.       loop set!(s*elt!) end end;      
  115.    
  116.    scale_col_by(c: INT, s: FLT) is
  117.       -- Scale the "c"th column of self by s
  118.       loop set_col!(c,s*col_elt!(c)) end end;
  119.       
  120.    scale_row_by(r: INT, s: FLT) is
  121.       -- Scale the "r"th row of self by s
  122.       loop set_row!(r,s*row_elt!(r)) end end;
  123.  
  124.    to_zero is
  125.       -- Set self to the zero matrix.
  126.       loop set!(0.0) end end;
  127.  
  128.    to_identity is
  129.       -- Set self to the identity matrix. Add extra zeroes if non-square.
  130.       to_zero; loop i::=nr.min(nc).times!; [i,i]:=1.0 end end;
  131.    
  132.    plus(m1:SAME):SAME
  133.       -- The sum of self and m1.  -> Self+m1
  134.       pre same_shape_p(m1) is
  135.       res::=create(nr,nc); res.to_sum_of(self,m1); return(res); end;
  136.  
  137.    add(m1:SAME)
  138.       -- Set self to its sum with m1.  Self <- Self + m1
  139.       pre same_shape_p(m1) is
  140.       loop set!(elt!+m1.elt!) end end;
  141.    
  142.    add_scalar(s: FLT) is
  143.       -- Add scalar "s" to all elements fo the matrix
  144.       loop set!(elt!+s) end;  end;
  145.    
  146.    add_scaled(scale:FLT, m1: SAME)
  147.       -- Set self to self + scale*m1
  148.       pre same_shape_p(m1) is
  149.       loop set!(elt!+scale*m1.elt!) end end;
  150.    
  151.    scale_and_add_scaled(scale_self,scale_m1:FLT, m1: SAME)
  152.       -- Set self to self_scale*self + scale_m1*m1
  153.       pre same_shape_p(m1) is
  154.       loop set!(scale_self*elt!+scale_m1*m1.elt!) end end;
  155.    
  156.    to_sum_of(m1,m2:SAME)
  157.       -- Set self to sum of m1 and m2. Self <- m1+m2
  158.       pre same_shape_p(m1) and same_shape_p(m2) is
  159.       loop set!(m1.elt!+m2.elt!) end end;
  160.    
  161.    minus(m1:SAME):SAME
  162.       -- ->Self-m1.
  163.       pre same_shape_p(m1) is
  164.       res::=create(nr,nc); loop res.set!(elt!-m1.elt!) end;
  165.       return(res);
  166.       end;
  167.    
  168.    to_difference_with(m1:SAME)
  169.       -- Self <- Self - m1.
  170.       pre same_shape_p(m1) is
  171.       loop set!(elt!-m1.elt!) end end;
  172.  
  173.    to_difference_of(m1,m2:SAME)
  174.       -- Set self to difference of m1 and m2. Self <- m1-m2
  175.       pre same_shape_p(m1) and same_shape_p(m2) is
  176.       loop set!(m1.elt!-m2.elt!) end end;
  177.  
  178.    add_difference_of(m1,m2: SAME)
  179.       -- Add to self the difference of m1 and m2. Self <- Self + m1-m2
  180.       pre same_shape_p(m1) and same_shape_p(m2) is
  181.       loop set!(elt! + (m1.elt!-m2.elt!)) end end;
  182.    
  183.    times(m1:SAME):SAME
  184.       -- Self times m1. -> Self*m1
  185.       pre m1.nr=nc is
  186.       res::=create(nr,m1.nc);
  187.       res.to_product_of(self,m1);
  188.       return(res);
  189.       end;
  190.    
  191.    to_product_of(a,b:SAME)
  192.       -- Set self to the matrix product of `a' and `b'.
  193.       pre nr=a.nr and nc=b.nc and a.nc=b.nr is
  194.       loop t::=inds!; x::=0.0;
  195.      a_row ::= t.t1; b_col ::= t.t2;
  196.      loop x:=x+a.row_elt!(a_row)*b.col_elt!(b_col) end;
  197.      [a_row,b_col]:= x end end;
  198.  
  199.    add_product_of(a,b: SAME)
  200.       -- Add the matrix product of `a' and `b' to self
  201.          -- Self <- Self+a*b
  202.       pre nr=a.nr and nc=b.nc and a.nc=b.nr is
  203.       loop t::=inds!; x::=0.0;
  204.      a_row ::= t.t1; b_col ::= t.t2;
  205.      loop x:=x+a.row_elt!(a_row)*b.col_elt!(b_col) end;
  206.      [a_row,b_col]:= [a_row,b_col]+x end end;
  207.    
  208.  
  209.    scale_and_add_product_of(scale_self,scale_prod: FLT, a,b: SAME)
  210.          -- Add the matrix product of `a' and `b' to scaled*self
  211.          -- SELF <- scale_self*SELF +scale_product*A*B
  212.      -- Blas 3 routine: might want to write other routines to call this
  213.       pre nr=a.nr and nc=b.nc and a.nc=b.nr is
  214.       loop t::=inds!; x::=0.0;
  215.      a_row ::= t.t1; b_col ::= t.t2;
  216.      loop x:=x+a.row_elt!(a_row)*b.col_elt!(b_col) end;
  217.      [a_row,b_col]:= scale_prod*[a_row,b_col]+scale_self*x end end;
  218.    
  219.    to_product_with_diagonal(v1:VEC) is
  220.       -- Set self to be the product of itself with a diagonal matrix whose
  221.       -- diagonal entries are the components of v1 truncated or extended
  222.       -- with zeroes to be the correct size.
  223.       -- Scale the columns of self with the elements of v1
  224.       -- Self <- Self*v1 (as diagonal of a matrix)
  225.       loop c::=nc.min(v1.dim).times!;
  226.      loop set_col!(c,v1[c]*col_elt!(c)) end end;
  227.       loop c::=v1.dim.upto!(nc-1); loop set_col!(c,0.0) end end end;
  228.  
  229.    times(v1:VEC):VEC is
  230.       -- A new vector equal to self times v1. -> self*v1
  231.       res ::= #VEC(nr);  times_in(v1,res);   return(res) end;
  232.  
  233.    times_in(v,vr:VEC)
  234.       -- Set `vr' to the result of self acting on `v'. vr <- self*v1
  235.       pre v.dim=nc and vr.dim=nr is
  236.       loop r::=vr.ind!; x::=0.0; loop x:=x+row_elt!(r)*v.elt! end;
  237.      vr[r]:=x end end;
  238.    
  239.    times_accum_in(v,vr: VEC)
  240.       -- Add to  `vr' the result of self acting on `v'. vr <- vr+self*v1
  241.       pre v.dim=nc and vr.dim=nr is
  242.       loop r::=vr.ind!; x::=0.0; loop x:=x+row_elt!(r)*v.elt! end;
  243.      vr[r]:=v[r]+x end end;
  244.  
  245.    transpose_times(v:VEC):VEC is
  246.       -- The transpose of self times `v' .
  247.       -- Used by vec*matrix
  248.       res ::= #VEC(nc);
  249.       transpose_times_in(v,res);
  250.       return(res) end;
  251.    
  252.    transpose_times_in(v,vr:VEC)
  253.       -- Set `vr' to the result of the transpose of self acting on `v'.
  254.       pre v.dim=nr and vr.dim=nc is
  255.       vr.to_zero;
  256.       i ::=0; loop until!(i=nc);
  257.      j::=0; loop until!(j=nr); vr[i]:=vr[i]+[j,i]*v[j]; j:=j+1 end;
  258.      i:=i+1
  259.      end; --loop
  260.       end; 
  261.    
  262.    transpose_times_add_into(v,vr: VEC)
  263.       -- Add to `vr' to the result of the transpose of self acting on `v'.
  264.       pre v.dim=nr and vr.dim=nc is
  265.       i ::=0; loop until!(i=nc);
  266.      j::=0; loop until!(j=nr); vr[i]:=vr[i]+[j,i]*v[j]; j:=j+1 end;
  267.      i:=i+1
  268.      end; --loop
  269.       end; 
  270.    
  271.    col(j:INT):VEC is
  272.       -- The `j'th col of self as a vector
  273.       v::=#VEC(nr); col_in(j,v); return(v) end;
  274.  
  275.    col_in(j:INT, v:VEC) is
  276.       -- Copy as much of the `j'th col of self into `v' as will fit.
  277.       loop v.set!(col_elt!(j)) end end;
  278.  
  279.    set_col(j:INT, v:VEC) is
  280.       -- Copy as much of `v' as will fit into the `j'th col of self 
  281.       loop set_col!(j,v.elt!) end; end;
  282.    
  283.    row(i:INT):VEC is
  284.       -- The `i'th row of self.
  285.       v::=#VEC(nc);  row_in(i,v); return(v) end;
  286.  
  287.    row_in(i:INT, v:VEC) is
  288.       -- Copy as much of the `i'th row of self into `v' as will fit.
  289.       loop v.set!(row_elt!(i)) end; end;
  290.    
  291.    set_row(i:INT, v:VEC) is
  292.       -- Copy as much of `v' as will fit into the `i'th row of self 
  293.       loop set_row!(i,v.elt!) end end;
  294.  
  295.    diagonal_in(v:VEC)
  296.       -- Copy as much of the diagonal of self into `v' as will fit.
  297.       pre v.dim = nr.min(nc) is
  298.       loop v.set!(diag_elt!) end; end;
  299.  
  300.    trace:FLT is
  301.       -- The trace of self. Sum of diagonal even if not square.
  302.       res ::= 0.0; loop res := res + diag_elt! end end;
  303.  
  304.    to_uniform_random is 
  305.       -- Make self's entries uniform in `[0.,1.)' and return it;
  306.       loop r::=row_ind!; loop [r,col_ind!] := (RND::uniform).flt; end; end;
  307.       end; -- to_uniform_random
  308.  
  309.    to_outer_product_of(u,v:VEC)
  310.       -- Set `self[r,c]' to `u.r * v.c' (like `u v^T').
  311.       pre nr=u.dim and nc=v.dim is
  312.       loop r_ind ::= u.ind!;
  313.      loop c_ind ::= v.ind!;
  314.         [r_ind,c_ind] := u[r_ind]*v[c_ind];
  315.         end; end; end;
  316.    
  317.    add_outer_product_of(u,v:VEC)
  318.       -- Add to  `self[r,c]'  `u.r * v.c' (like `u v^T').
  319.       pre nr=u.dim and nc=v.dim is
  320.       loop r_ind ::= u.ind!;
  321.      loop c_ind ::= v.ind!;
  322.         [r_ind,c_ind] := [r_ind,c_ind] + u[r_ind]*v[c_ind];
  323.         end; end; end;
  324.  
  325.  
  326. end; -- class MAT
  327. --~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  328. class TEST_MAT is
  329.    -- Test of `MAT' 
  330.    include TEST;
  331.    
  332.    main is
  333.       -- Test of `MAT' and `NR_SVD'.
  334.       class_name("MAT");
  335.       m ::= #MAT(||1.1,4.2|,|3.7,5.6|,|5.9,6.3||);
  336.       test("nr", m.nr.str, "3");
  337.       test("nc", m.nc.str, "2");
  338.       test("create_from_s and to_s", m.str,
  339.          "||1.1,4.2|,|3.7,5.6|,|5.9,6.3||");
  340.         -- Having tested m.str, we should be able to
  341.         -- use it in the tests that follow to avoid
  342.         -- heavy dependance on the exact printout format of flts
  343.       m2 ::=#MAT(3,2);
  344.       m2.to(m);
  345.       test("to and create", m2.str, m.str);
  346.       m3 ::= #MAT(2,1);
  347.       
  348.       m3.to_portion_of(m);
  349.       m_portion_test ::= #MAT(||1.10|,|3.70||);
  350.       test("to_portion_of", m3.str, m_portion_test.str);
  351.       
  352.       m4 ::= m.copy;
  353.       test("copy",m4.str,m.str);
  354.       
  355.       m4.to_zero;
  356.       m_test ::= #MAT(||0.0,0.0|,|0.0,0.0|,|0.0,0.0||);
  357.       test("to_zero", m4.str,m_test.str);
  358.  
  359.       m_test.copy(||1.0,0.0|,|0.0,1.0|,|0.0,0.0||);
  360.       m4.to_identity;
  361.       test("to_identity", m4.str,m_test.str);
  362.       
  363.       m_test.copy(||2.10,4.20|,|3.70,6.60|,|5.90,6.30||);
  364.       test("plus", (m4+m).str, m_test.str);
  365.       
  366.       m_test.copy(||2.10, 4.20|, |3.70, 6.60|, |5.90, 6.30||);
  367.       m5 ::= m4.copy;
  368.       m5.add(m);
  369.       test("add", m5.str,m_test.str);
  370.  
  371.       m_test.copy(||0.10, 4.20|, |3.70, 4.60|, |5.90, 6.30||);
  372.       test("minus", (m-m4).str,m_test.str);
  373.       
  374.       m_test.copy(||2.0,0.0|,|0.0,2.0|,|0.0,0.0||);
  375.       m6 ::= m4.copy;
  376.       m6.scale_by(2.0);
  377.       test("scale_by", m6.str,m_test.str);
  378.       
  379.       m_test.copy(||1.10, 4.20|, |3.70, 5.60|, |5.90, 6.30||);
  380.       m_ident ::= MAT::create_identity(2);
  381.       test("times", (m*m_ident).str, m_test.str);
  382.  
  383.       v ::= #VEC(2);
  384.       v.to_ones;
  385.       m.to_product_with_diagonal(v);
  386.       test("times", m.str, m_test.str);
  387.  
  388.       m_trans_test ::= #MAT(||1.10, 3.70, 5.90|,|4.20,5.60,6.30||);
  389.       m_trans::=m.transpose;
  390.       test("transpose",m_trans.str,m_trans_test.str);
  391.       
  392.       m_trans.to_zero;
  393.       m_trans_test.copy(||0.0, 0.00, 0.00|,|0.0,0.0,0.0||);
  394.       test("to_zero",m_trans.str,m_trans_test.str);
  395.  
  396.       v.to_ones;
  397.       m_mul ::= m*v;
  398.       v_res ::= #VEC(|5.3,9.3,12.2|);
  399.       test("times vector",m_mul.str,v_res.str);
  400.  
  401.       v1 ::= #VEC(|5.0,7.0|);
  402.       v2 ::= #VEC(|3.0,9.0,4.0|);
  403.       m_res ::= v1.outer_product(v2);
  404.       m_trans.copy(||15.0,45.0,20.0|,|21.0,63.0,28.0||);
  405.       test("outer product (VEC test)", m_res.str,m_trans.str);
  406.       
  407.       m_trans := m.transpose;
  408.       v_m ::= v*m_trans;
  409.       v_test ::= #VEC(|5.3,9.3,12.2|);
  410.       test("vector * matrix i.e. transpose_times",v_m.str,v_test.str);
  411.       
  412. --      test("transpose", m.transpose.str,
  413. --     "((1.10, 3.70, 5.90), (4.20, 5.60, 6.30))");
  414. --      test("to_transpose_of", MAT::create(2,3).to_transpose_of(m).str,
  415. --     "((1.10, 3.70, 5.90), (4.20, 5.60, 6.30))");
  416. --      test("transpose_act_on", m.transpose.transpose_act_on(v).str,
  417. --     "(1.10, 3.70, 5.90)");
  418. --      m.transpose.transpose_act_on_in(v,v2);
  419. --      test("transpose_act_on_in", v2.str, "(1.10, 3.70, 5.90)");
  420.  
  421. --      v:VEC :=VEC::create_from_s("(1.0, 0.0)");
  422. --      test("act_on", m.act_on(v).str, "(1.10, 3.70, 5.90)");
  423. --      v2:VEC :=VEC::create(3); m.act_on_in(v,v2);
  424. --      test("act_on_in", v2.str, "(1.10, 3.70, 5.90)");
  425. --      v3:VEC :=VEC::create_from_s("(0.0)");
  426. --      test("column", m.column(0).str, "(1.10, 3.70, 5.90)");
  427. --      m.column_in(0,v2);
  428. --      test("column_in", v2.str, "(1.10, 3.70, 5.90)");
  429. --      test("column_from", m.copy.column_from(0,v2.to_zero).str,
  430. --     "((0.00, 4.20), (0.00, 5.60), (0.00, 6.30))");
  431. --      test("row", m.row(0).str, "(1.10, 4.20)");
  432. --      m.row_in(0,v);
  433. --      test("row_in", v.str, "(1.10, 4.20)");
  434. --      test("row_from", m.copy.row_from(0,v.to_zero).str,
  435. --     "((0.00, 0.00), (3.70, 5.60), (5.90, 6.30))");
  436. --      m.diagonal_in(v);
  437. --      test("diagonal_in", v.str, "(1.10, 5.60)");
  438. --      test("trace", m.trace.str, (m[0,0]+m[1,1]).str);
  439. --      a:MAT:=MAT::create(3,4).to_uniform_random;
  440. --      svd_mat_test(1, a);
  441. --      a:=MAT::create(4,3).to_uniform_random;
  442. --      a:=MAT::create_from_s("((0.,1.,0.),(0.,1.,1.),(0.,0.,0.))");
  443. --      test("to_uniform_random",
  444. --     lfm.copy.to_uniform_random.nr.str, "3");
  445. --      test("to_outer_product_of",
  446. --     em.to_outer_product_of
  447. --           (VEC::create_from_s("(1.,1.,1.)"),
  448. --          VEC::create_from_s("(1.,1.,1.)")).str,
  449. --     "((1.00, 1.00, 1.00), (1.00, 1.00, 1.00), (1.00, 1.00, 1.00))");
  450.       finish;
  451.    end; -- main
  452.  
  453.  
  454. end; -- class TEST_MAT
  455.  
  456. --~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  457.