Mercurial > hg > octave-max
changeset 4:b4df021f796c
[project @ 1993-08-08 01:26:08 by jwe]
Initial revision
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/scripts/Makefile.in @@ -0,0 +1,80 @@ +# +# Makefile for octave's scripts directory +# +# John W. Eaton +# jwe@che.utexas.edu +# Department of Chemical Engineering +# The University of Texas at Austin + +TOPDIR = .. + +srcdir = @srcdir@ +VPATH = @srcdir@ + +include $(TOPDIR)/Makeconf + +SOURCES = computer.in *.m + +DISTFILES = Makefile.in $(SOURCES) + +MFILES = $(wildcard $(srcdir)/*.m) +MFILES_NO_DIR = $(notdir $(MFILES)) + +all: computer.m +.PHONY: all + +computer.m: computer.in + if test "$(target_host_type)" = unknown ; then \ + sed -e "s/%target_host_type%/Hi Dave, I'm a HAL-9000/" \ + $(srcdir)/computer.in > computer.m ; \ + else \ + sed -e "s/%target_host_type%/$(target_host_type)/" \ + $(srcdir)/computer.in > computer.m ; \ + fi + +check: all +.PHONY: check + +install: all + if test -d $(libsubdir) ; then true ; \ + else $(TOPDIR)/mkpath $(libsubdir) ; fi + for f in $(MFILES_NO_DIR) ; do \ + rm -f $(libdir)/$$f ; \ + $(INSTALL_DATA) $(srcdir)/$$f $(libsubdir)/$$f ; \ + done + rm -f $(libsubdir)/computer.m + $(INSTALL_DATA) computer.m $(libsubdir)/computer.m +.PHONY: install + +uninstall: + for f in $(MFILES_NO_DIR) ; do rm -f $(libsubdir)/$$f ; done +.PHONY: uninstall + +clean: + rm -f computer.m +.PHONY: clean + +tags: $(SOURCES) + ctags $(SOURCES) + +TAGS: $(SOURCES) + etags $(SOURCES) + +mostlyclean: clean +.PHONY: mostlyclean + +distclean: clean + rm -f Makefile +.PHONY: distclean + +realclean: distclean + rm -f tags TAGS +.PHONY: realclean + +local-dist: + ln $(DISTFILES) ../`cat ../.fname`/scripts +.PHONY: local-dist + +dist: + ln $(DISTFILES) ../`cat ../.fname`/scripts +.PHONY: dist
new file mode 100644 --- /dev/null +++ b/scripts/general/columns.m @@ -0,0 +1,15 @@ +function nc = columns (x) + +# usage: columns (x) +# +# Return the the number of columns in x. +# +# See also: size, rows, length, is_scalar, is_vector, is_matrix + + if (nargin != 1) + error ("usage: columns (x)"); + endif + + [nr, nc] = size (x); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/fliplr.m @@ -0,0 +1,17 @@ +function y = fliplr (x) + +# usage: fliplr (x) +# +# Return x with the columns swapped. +# +# See also: flipu, rot90 + + if (nargin != 1) + error ("usage: fliplr (x)"); + endif + + y = x; + nc = columns (x); + y = x (:, nc:-1:1); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/flipud.m @@ -0,0 +1,17 @@ +function y = flipud (x) + +# usage: flipud (x) +# +# Return x with the rows swapped. +# +# See also: fliplr, rot90 + + if (nargin != 1) + error ("usage: flipud (x)"); + endif + + y = x; + nr = rows (x); + y = x (nr:-1:1, :); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/int2str.m @@ -0,0 +1,19 @@ +function retval = int2str (x) + +# usage: int2str (x) +# +# Round x to the nearest integer and format as a string. +# +# See also: sprintf, num2str + + if (nargin == 1) + if (rows (x) == 1 && columns (x) == 1) + retval = sprintf ("%f\n", round (x)); + else + error ("int2str: expecting scalar argument"); + endif + else + error ("usage: int2str (x)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/is_matrix.m @@ -0,0 +1,17 @@ +function retval = is_matrix (x) + +# usage: is_matrix (x) +# +# Return 1 if the number of rows and columns of x are both greater +# than 1. +# +# See also: size, rows, columns, length, is_scalar, is_vector + + if (nargin == 1) + [nr, nc] = size (x); + retval = (nr > 1 && nc > 1); + else + error ("usage: is_matrix (x)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/is_scalar.m @@ -0,0 +1,16 @@ +function retval = is_scalar (x) + +# usage: is_scalar (x) +# +# Return 1 if the number of rows and columns of x are both equal to 1. +# +# See also: size, rows, columns, length, is_scalar, is_matrix + + if (nargin == 1) + [nr, nc] = size (x); + retval = (nr == 1 && nc == 1); + else + error ("usage: is_scalar (x)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/is_vector.m @@ -0,0 +1,17 @@ +function retval = is_vector (x) + +# usage: is_vector (x) +# +# Return 1 if the either the number of rows (columns) of x is 1 and +# the number of columns (rows) is greater than one. Otherwise, return 0. +# +# See also: size, rows, columns, length, is_scalar, is_matrix + + if (nargin == 1) + [nr, nc] = size (x); + retval = ((nr == 1 && nc > 1) || (nc == 1 && nr > 1)); + else + error ("usage: is_vector (x)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/isempty.m @@ -0,0 +1,13 @@ +function retval = isempty (var) + +# usage: isempty (x) +# +# Return 1 if the argument is an empty matrix. Otherwise, return 0. + + if (nargin != 1) + error ("usage: isempty (var)"); + endif + + retval = (rows (var) == 0 || columns (var) == 0); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/length.m @@ -0,0 +1,15 @@ +function len = length (x) + +# usage: length (x) +# +# Return the number of rows or columns, whichever is greater. +# +# See also: size, rows, columns, is_scalar, is_vector, is_matrix + + if (nargin != 1) + error ("usage: length (x)"); + endif + + len = max (size (x)); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/linspace.m @@ -0,0 +1,40 @@ +function retval = linspace (x1, x2, n) + +# usage: linspace (x1, x2, n) +# +# Return a vector of n equally spaced points between x1 and x2 +# inclusive. +# +# If the final argument is omitted, n = 100 is assumed. +# +# All three arguments must be scalars. +# +# See also: logspace + + if (nargin == 2) + npoints = 100; + elseif (nargin == 3) + if (length (n) == 1) + npoints = n; + else + error ("linspace: arguments must be scalars"); + endif + else + error ("usage: linspace (x1, x2 [, n])"); + endif + + if (npoints < 2) + error ("linspace: npoints must be greater than 2"); + endif + + if (length (x1) == 1 && length (x2) == 1) + delta = (x2 - x1) / (npoints - 1); + retval = zeros (1, npoints); + for i = 0:npoints-1 + retval (i+1) = x1 + i * delta; + endfor + else + error ("linspace: arguments must be scalars"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/logspace.m @@ -0,0 +1,50 @@ +function retval = logspace (x1, x2, n) + +# usage: logspace (x1, x2, n) +# +# Return a vector of n logarithmically equally spaced points between +# x1 and x2 inclusive. +# +# If the final argument is omitted, n = 50 is assumed. +# +# All three arguments must be scalars. +# +# Note that if if x2 is pi, the points are between 10^x1 and pi, NOT +# 10^x1 and 10^pi. +# +# Yes, this is pretty stupid, because you could achieve the same +# result with logspace (x1, log10 (pi)), but Matlab does this, and +# claims that is useful for signal processing applications. +# +# See also: linspace + + if (nargin == 2) + npoints = 50; + elseif (nargin == 3) + if (length (n) == 1) + npoints = n; + else + error ("logspace: arguments must be scalars"); + endif + else + error ("usage: logspace (x1, x2 [, n])"); + endif + + if (npoints < 2) + error ("logspace: npoints must be greater than 2"); + endif + + if (length (x1) == 1 && length (x2) == 1) + x2_tmp = x2; + if (x2 == pi) + x2_tmp = log10 (pi); + endif + retval = linspace (x1, x2_tmp, npoints); + for i = 1:npoints + retval(i) = 10 ^ retval(i); + endfor + else + error ("logspace: arguments must be scalars"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/num2str.m @@ -0,0 +1,19 @@ +function retval = num2str (x) + +# usage: num2str (x) +# +# Format x as a string. +# +# See also: sprintf, int2str + + if (nargin == 1) + if (rows (x) == 1 && columns (x) == 1) + retval = sprintf ("%g", x); + else + error ("num2str: expecting scalar argument"); + endif + else + error ("usage: num2str (x)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/perror.m @@ -0,0 +1,92 @@ +function perror (name, err) + +# usage: perror (name, err) +# +# Print an error message for error number `err' from function "name". +# +# Messages correspond to the following subroutine versions: +# +# npsol : 4.0 +# qpsol : 3.2 + + if (nargin != 2) + error ("usage: perror (name, err)"); + endif + + if (! isstr (name)) + error ("perror: first argument must be a string"); + endif + + if (! is_scalar (err)) + error ("perror: second argument must be a scalar"); + endif + + if (strcmp (name, "fsolve")) + + if (info == -1) + printf ("input error\n"); + elseif (info == 1) + printf ("solution converged to requested tolerance\n"); + elseif (info == 4) + printf ("iteration limit exceeded\n"); + elseif (info == 3) + printf ("iteration is not making good progress\n"); + else + error ("perror: unrecognized error code for fsolve"); + endif + + elseif (strcmp (name, "npsol")) + + if (err == 0) + printf ("optimal solution found\n"); + elseif (err == 1) + printf ("weak local solution found\n"); + elseif (err == 2) + printf ("no feasible point for linear constraints and bounds\n"); + elseif (err == 3) + printf ("no feasible point found for nonlinear constraints\n"); + elseif (err == 4) + printf ("iteration limit reached\n"); + elseif (err == 6) + printf ("current point cannot be improved upon\n"); + elseif (err == 7) + printf ("user-supplied derivatives appear to be incorrect\n"); + elseif (err == 9) + printf ("internal error: invalid input parameter\n"); + else + error ("perror: unrecognized error code for npsol"); + endif + + elseif (strcmp (name, "qpsol")) + + if (err == 0) + printf ("optimal solution found\n"); + elseif (err == 1) + printf ("weak local solution found\n"); + elseif (err == 2) + printf ("solution appears to be unbounded\n"); + elseif (err == 3) + printf ("solution appears optimal, but optimality can't be verified\n"); + elseif (err == 4) + printf ("iterates of the QP phase appear to be cycling\n"); + elseif (err == 5) + printf ("iteration limit reached during QP phase\n"); + elseif (err == 6) + printf ("no feasible point found during LP phase\n"); + elseif (err == 7) + printf ("iterates of the LP phase appear to be cycling\n"); + elseif (err == 8) + printf ("iteration limit reached during LP phase\n"); + elseif (err == 9) + printf ("internal error: invalid input parameter\n"); + else + error ("perror: unrecognized error code for qpsol"); + endif + + else + + error ("perror: unrecognized function name"); + + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/rem.m @@ -0,0 +1,28 @@ +function retval = rem (x, y) + +# usage: rem (x, y) +# +# Return remainder (x, y). + + if (nargin != 2) + error ("usage: rem (x, y)"); + endif + + if (any (size (x) != size (y))) + error ("rem: argument sizes must agree") + endif + +# Matlab allows complex arguments, but as far as I can tell, that's a +# bunch of hooey. + + if (any (any (imag (x))) || any (any (imag (y)))) + error ("rem: complex arguments are not allowed"); + endif + + if (nargin == 2) + retval = x - y .* fix (x ./ y); + else + error ("usage: rem (x, y)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/reshape.m @@ -0,0 +1,25 @@ +function retval = reshape (a, m, n) + +# usage: reshape (a, m, n) +# +# Form an m x n matrix from the elements of a (taken in Fortran's +# column major ordering). +# +# See also: `:', do_fortran_indexing + + if (nargin != 3) + error ("usage: reshape (a, m, n)"); + else + [nr, nc] = size (a); + if (nr * nc == m * n) + tmp = do_fortran_indexing; + do_fortran_indexing = "true"; + retval = zeros (m, n); + retval (:) = a; + do_fortran_indexing = tmp; + else + error ("reshape: sizes must match"); + endif + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/rot90.m @@ -0,0 +1,39 @@ +function y = rot90 (x, k) + +# usage: rot90 (x, k) +# +# Rotate the matrix x counterclockwise k*90 degrees. +# +# If the second argument is omitted, k is taken to be 1. +# +# See also: flipud, fliplr + + if (nargin < 2) + k = 1; + endif + + if (imag (k) != 0 || fix (k) != k) + error ("rot90: k must be an integer"); + endif + + if (nargin == 1 || nargin == 2) + k = rem (k, 4); + if (k < 0) + k = k + 4; + endif + if (k == 0) + y = x; + elseif (k == 1) + y = flipud (x'); + elseif (k == 2) + y = flipud (fliplr (x)); + elseif (k == 3) + y = (flipud (x))'; + else + error ("rot90: internal error!"); + endif + else + error ("usage: rot90 (x [, k])"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/rows.m @@ -0,0 +1,15 @@ +function nr = rows (x) + +# usage: rows (x) +# +# Return the the number of rows in x. +# +# See also: size, columns, length, is_scalar, is_vector, is_matrix + + if (nargin != 1) + error ("usage: rows (x)"); + endif + + [nr, nc] = size (x); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/tril.m @@ -0,0 +1,32 @@ +function retval = tril (x, k) + +# usage: triu (x, k) +# +# Return the lower triangular part of x above the k-th diagonal. If +# the second argument is omitted, k = 0 is assumed. +# +# See also: triu, diag + + if (nargin > 0) + [nr, nc] = size (x); + retval = x; + endif + + if (nargin == 1) + k = 0; + elseif (nargin == 2) + max_nr_nc = max (nr, nc); + if ((k > 0 && k > nr - 1) || (k < 0 && k < 1 - nc)) + error ("tril: requested diagonal out of range") + endif + else + error ("usage: tril (x [, k])"); + endif + + for i = 1:nr + for j = i+1-k:nc + retval (i, j) = 0.0; + endfor + endfor + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/general/triu.m @@ -0,0 +1,32 @@ +function retval = triu (x, k) + +# usage: triu (x, k) +# +# Return the upper triangular part of x above the k-th diagonal. If +# the second argument is omitted, k = 0 is assumed. +# +# See also: tril, diag + + if (nargin > 0) + [nr, nc] = size (x); + retval = x; + endif + + if (nargin == 1) + k = 0; + elseif (nargin == 2) + max_nr_nc = max (nr, nc); + if ((k > 0 && k > nc - 1) || (k < 0 && k < 1 - nr)) + error ("triu: requested diagonal out of range") + endif + else + error ("usage: triu (x [, k])"); + endif + + for j = 1:nc + for i = j+1-k:nr + retval (i, j) = 0.0; + endfor + endfor + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/cond.m @@ -0,0 +1,27 @@ +function retval = cond (a) + +# usage: cond (a) +# +# Return the condition number of a, computed using the singular values +# of a. +# +# See also: norm, svd + + if (nargin == 1) + [nr, nc] = size (a); + if (nr == 0 && nc == 0) + if (strcmp (propagate_empty_matrices, "false")) + error ("cond: empty matrix is invalid as argument"); + endif + if (strcmp (propagate_empty_matrices, "warn")) + printf ("warning: cond: argument is empty matrix\n"); + endif + retval = 0.0; + endif + sigma = svd (a); + retval = sigma (1) / sigma (length (sigma)); + else + error ("usage: cond (a)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/norm.m @@ -0,0 +1,83 @@ +function retval = norm (x, p) + +# usage: norm (x, p) +# +# Compute the p-norm of x. +# +# If x is a matrix: +# +# value of p norm returns +# ---------- ------------ +# 1 1-norm, the largest column sum of x +# 2 largest singular value of x +# Inf infinity norm, the largest row sum of x +# "fro" Frobenius norm of x, sqrt (sum (diag (x' * x))) +# +# If x is a vector or a scalar: +# +# value of p norm returns +# ---------- ------------ +# Inf max (abs (x)) +# -Inf min (abs (x)) +# other p-norm of x, sum (abs (x) .^ p) ^ (1/p) +# +# If the second argument is missing, p = 2 is assumed. +# +# See also: cond, svd + + if (nargin < 1 || nargin > 2) + error ("usage: norm (x [, p])") + endif + +# Do we have a vector or matrix as the first argument? + + if (rows (x) == 1 || columns (x) == 1) + + if (nargin == 2) + if (isstr (p)) + if (strcmp (p, "fro")) + retval = sqrt (sum (diag (x' * x))); + else + error ("norm: unrecognized norm"); + endif + else + if (p == Inf) + retval = max (abs (x)); + elseif (p == -Inf) + retval = min (abs (x)); + else + retval = sum (abs (x) .^ p) ^ (1/p); + endif + endif + elseif (nargin == 1) + retval = sum (abs (x) .^ 2) ^ 0.5; + endif + + else + + if (nargin == 2) + if (isstr (p)) + if (strcmp (p, "fro")) + retval = sqrt (sum (diag (x' * x))); + else + error ("norm: unrecognized norm"); + endif + else + if (p == 1) + retval = max (sum (abs (real (x)) + abs (imag (x)))); + elseif (p == 2) + s = svd (x); + retval = s (1); + elseif (p == Inf) + xp = x'; + retval = max (sum (abs (real (xp)) + abs (imag (xp)))); + endif + endif + elseif (nargin == 1) + s = svd (x); + retval = s (1); + endif + + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/rank.m @@ -0,0 +1,25 @@ +function retval = rank (A, tol) + +# usage: rand (a, tol) +# +# Return the rank of the matrix a. The rank is taken to be the number +# of singular values of a that are greater than tol. +# +# If the second argument is omitted, it is taken to be +# +# tol = max (size (a)) * sigma (1) * eps; +# +# where eps is machine precision and sigma is the largest singular +# value of a. + + if (nargin == 1) + sigma = svd (A); + tolerance = max (size (A)) * sigma (1) * eps; + elseif (nargin == 2) + tolerance = tol; + else + error ("usage: rank (A)"); + endif + retval = sum (sigma > tolerance); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/trace.m @@ -0,0 +1,13 @@ +function y = trace (x) + +# usage: trace (x) +# +# Returns the trace (the sum of the diagonal elements) of x. + + if (nargin != 1) + error ("usage: trace (x)"); + endif + + y = sum (diag (x)); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/miscellaneous/menu.m @@ -0,0 +1,47 @@ +function s = menu (t,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16) + +# usage: menu (title, opt1, opt2, ..., opt16) +# +# See also: disp, printf, input + + if (nargin < 2 || nargin > 17) + error ("usage: menu (title, opt1, opt2, ..., opt16)"); + endif + +# Force pending output to appear before the menu. + + fflush (stdout); + +# Don't send the menu through the pager since doing that can cause +# major confusion. + + save_page_screen_output = page_screen_output; + page_screen_output = "false"; + + if (! isempty (t)) + disp (t); + printf ("\n"); + endif + + nopt = nargin - 1; + + s = 0; + while (1) + page_screen_output = "false"; + for i = 1:nopt + command = sprintf ("printf (\" [%2d] \"); disp (x%d)", i, i); + eval (command); + endfor + printf ("\n"); + page_screen_output = save_page_screen_output; + s = input ("pick a number, any number: "); + if (s < 1 || s > nopt) + printf ("\nerror: input out of range\n\n"); + else + break; + endif + endwhile + + page_screen_output = save_page_screen_output; + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/miscellaneous/texas_lotto.m @@ -0,0 +1,34 @@ +function picks = texas_lotto () + +# usage: texas_lotto +# +# Pick 6 unique numbers between 1 and 50 that are guaranteed to win +# the Texas Lotto. +# +# See also: rand + + if (nargin != 0) + disp ("win_texas_lotto: ignoring extra arguments"); + endif + + picks = zeros (1,6); + picks (1) = round (50-49*(1-rand)); + n = 2; + while (n < 7) + tmp = round (50-49*(1-rand)); + equal = 0; + for i = 1:n + if (tmp == picks (i)) + equal = 1; + break; + endif + endfor + if (! equal) + picks (n) = tmp; + n++; + endif + endwhile + + picks = sort (picks); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/__plr__.m @@ -0,0 +1,108 @@ +function polar_int (theta, rho) + + if (nargin == 1) + [nr, nc] = size (theta); + if (nr == 1) + theta = theta'; + tmp = nr; + nr = nc; + nc = tmp; + endif + theta_i = imag (theta); + if (theta_i) + rho = theta_i; + theta = real (theta); + else + rho = theta; + theta = (1:nr)'; + endif + endif + + if (nargin <= 2) + if (imag (theta)) + theta = real (theta); + endif + if (imag (rho)) + rho = real (rho); + endif + if (is_scalar (theta)) + if (is_scalar (rho)) + x = rho * cos (theta); + y = rho * sin (theta); + plot_2_s_s (x, y); + endif + elseif (is_vector (theta)) + if (is_vector (rho)) + if (length (theta) != length (rho)) + error ("error: polar: vector lengths must match"); + endif + if (rows (rho) == 1) + rho = rho'; + endif + if (rows (theta) == 1) + theta = theta'; + endif + x = rho .* cos (theta); + y = rho .* sin (theta); + plot_2_v_v (x, y); + elseif (is_matrix (rho)) + [t_nr, t_nc] = size (theta); + if (t_nr == 1) + theta = theta'; + tmp = t_nr; + t_nr = t_nc; + t_nc = tmp; + endif + [r_nr, r_nc] = size (rho); + if (t_nr != r_nr) + rho = rho' + tmp = r_nr; + r_nr = r_nc; + r_nc = tmp; + endif + if (t_nr != r_nr) + error ("error: polar: vector and matrix sizes must match"); + endif + x = diag (cos (theta)) * rho; + y = diag (sin (theta)) * rho; + plot_2_v_m (x, y); + endif + elseif (is_matrix (theta)) + if (is_vector (rho)) + [r_nr, r_nc] = size (rho); + if (r_nr == 1) + rho = rho'; + tmp = r_nr; + r_nr = r_nc; + r_nc = tmp; + endif + [t_nr, t_nc] = size (theta); + if (r_nr != t_nr) + theta = rho' + tmp = t_nr; + t_nr = t_nc; + t_nc = tmp; + endif + if (r_nr != t_nr) + error ("error: polar: vector and matrix sizes must match"); + endif + diag_r = diag (r); + x = diag_r * cos (theta); + y = diag_r * sin (theta); + plot_2_m_v (x, y); + elseif (is_matrix (rho)) + if (size (rho) != size (theta)) + error ("error: polar: matrix dimensions must match"); + endif + x = rho .* cos (theta); + y = rho .* sin (theta); + plot_2_m_m (x, y); + endif + endif + else + usage = sprintf ("usage: polar_int (x)\n"); + usage = sprintf ("%s polar_int (x, y)", usage); + error (usage); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/__plt2mm__.m @@ -0,0 +1,25 @@ +function plot_2_m_m (x, y) + + if (nargin != 2) + error ("usage: plot_2_m_m (x, y)"); + endif + + [x_nr, x_nc] = size (x); + [y_nr, y_nc] = size (y); + + if (x_nr == y_nr && x_nc == y_nc) + if (x_nc > 0) + tmp = [x, y]; + command = sprintf ("gplot tmp(:,%d:%d:%d)", 1, x_nc, x_nc+1); + for i = 2:x_nc + command = sprintf ("%s, tmp(:,%d:%d:%d)", command, i, x_nc, x_nc+i); + endfor + eval (command); + else + error ("plot_2_m_m: arguments must be a matrices"); + endif + else + error ("plot_2_m_m: matrix dimensions must match"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/__plt2mv__.m @@ -0,0 +1,39 @@ +function plot_2_m_v (x, y) + + if (nargin != 2) + error ("usage: plot_2_m_v (x, y)"); + endif + + [x_nr, x_nc] = size (x); + [y_nr, y_nc] = size (y); + + if (y_nr == 1) + y = y'; + tmp = y_nr; + y_nr = y_nc; + y_nc = tmp; + endif + + if (x_nr == y_nr) + 1; + elseif (x_nc == y_nr) + x = x'; + tmp = x_nr; + x_nr = x_nc; + x_nc = tmp; + else + error ("plot_2_m_v: matrix dimensions must match"); + endif + + if (x_nc > 0) + tmp = [x, y]; + command = sprintf ("gplot tmp(:,%d:%d:%d)", 1, x_nc, x_nc+1); + for i = 2:x_nc + command = sprintf ("%s, tmp(:,%d:%d:%d)", command, i, x_nc-i+1, x_nc+1); + endfor + eval (command); + else + error ("plot_2_m_v: arguments must be a matrices"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/__plt2ss__.m @@ -0,0 +1,18 @@ +function plot_2_s_s (x, y) + + if (nargin != 2) + error ("usage: plot_2_s_s (x, y)"); + endif + + [x_nr, x_nc] = size (x); + [y_nr, y_nc] = size (y); + + if (x_nr == 1 && x_nr == y_nr && x_nc == 1 && x_nc == y_nc) + tmp = [x, y]; + command = sprintf ("gplot tmp"); + eval ("gplot tmp"); + else + error ("plot_2_s_s: arguments must be scalars"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/__plt2vm__.m @@ -0,0 +1,39 @@ +function plot_2_v_m (x, y) + + if (nargin != 2) + error ("usage: plot_2_v_m (x, y)"); + endif + + [x_nr, x_nc] = size (x); + [y_nr, y_nc] = size (y); + + if (x_nr == 1) + x = x'; + tmp = x_nr; + x_nr = x_nc; + x_nc = tmp; + endif + + if (x_nr == y_nr) + 1; + elseif (x_nr == y_nc) + y = y'; + tmp = y_nr; + y_nr = y_nc; + y_nc = tmp; + else + error ("plot_2_v_m: matrix dimensions must match"); + endif + + if (y_nc > 0) + tmp = [x, y]; + command = sprintf ("gplot tmp(:,%d:%d:%d)", 1, x_nc, x_nc+1); + for i = 2:y_nc + command = sprintf ("%s, tmp(:,%d:%d:%d)", command, 1, i, i+1); + endfor + eval (command); + else + error ("plot_2_v_m: arguments must be a matrices"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/__plt2vv__.m @@ -0,0 +1,31 @@ +function plot_2_v_v (x, y) + + if (nargin != 2) + error ("usage: plot_2_m_m (x, y)"); + endif + + [x_nr, x_nc] = size (x); + [y_nr, y_nc] = size (y); + + if (x_nr == 1) + x = x'; + tmp = x_nr; + x_nr = x_nc; + x_nc = tmp; + endif + + if (y_nr == 1) + y = y'; + tmp = y_nr; + y_nr = y_nc; + y_nc = tmp; + endif + + if (x_nr != y_nr) + error ("plot_2_v_v: vector lengths must match"); + endif + + tmp = [x, y]; + eval ("gplot tmp"); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/__plt__.m @@ -0,0 +1,51 @@ +function plot_int (x1, x2) + + if (nargin == 1) + [nr, nc] = size (x1); + if (nr == 1) + x1 = x1'; + tmp = nr; + nr = nc; + nc = tmp; + endif + x1_i = imag (x1); + if (x1_i) + x2 = x1_i; + x1 = real (x1); + else + x2 = x1; + x1 = (1:nr)'; + endif + endif + + if (nargin <= 2) + if (imag (x1)) + x1 = real (x1); + endif + if (imag (x2)) + x2 = real (x2); + endif + if (is_scalar (x1)) + if (is_scalar (x2)) + plot_2_s_s (x1, x2); + endif + elseif (is_vector (x1)) + if (is_vector (x2)) + plot_2_v_v (x1, x2); + elseif (is_matrix (x2)) + plot_2_v_m (x1, x2); + endif + elseif (is_matrix (x1)) + if (is_vector (x2)) + plot_2_m_v (x1, x2); + elseif (is_matrix (x2)) + plot_2_m_m (x1, x2); + endif + endif + else + usage = sprintf ("usage: plot_int (x)\n"); + usage = sprintf ("%s plot_int (x, y)", usage); + error (usage); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/bar.m @@ -0,0 +1,84 @@ +function [xb, yb] = bar (x, y) + +# usage: [xb, yb] = bar (x, y) +# +# Given two vectors of x-y data, bar produces a bar graph. +# +# If only one argument is given, it is taken as a vector of y-values +# and the x coordinates are taken to be the indices of the elements. +# +# If two output arguments are specified, the data are generated but +# not plotted. For example, +# +# bar (x, y); +# +# and +# +# [xb, yb] = bar (x, y); +# plot (xb, yb); +# +# are equivalent. +# +# See also: plot, semilogx, semilogy, loglog, polar, mesh, contour, +# stairs, gplot, gsplot, replot, xlabel, ylabel, title + + if (nargin == 1) + if (is_vector (x)) + len = 3 * length (x) + 1; + xb = yb = zeros (len, 1); + xb(1) = 0.5; + yb(1) = 0; + k = 1; + for i = 2:3:len + xb(i) = k-0.5; + xb(i+1) = k+0.5; + xb(i+2) = k+0.5; + yb(i) = x(k); + yb(i+1) = x(k); + yb(i+2) = 0.0; + k++; + endfor + else + error ("bar: argument must be a vector"); + endif + elseif (nargin == 2) + if (is_vector (x) && is_vector (y)) + xlen = length (x); + ylen = length (y); + if (xlen == ylen) + len = 3 * xlen + 1; + xb = yb = zeros (len, 1); + delta = (x(2) - x(1)) / 2.0; + xb(1) = x(1) - delta; + yb(1) = 0.0; + k = 1; + for i = 2:3:len + xb(i) = xb(i-1); + xb(i+1) = xb(i) + 2.0 * delta; + xb(i+2) = xb(i+1); + yb(i) = y(k); + yb(i+1) = y(k); + yb(i+2) = 0.0; + if (k < xlen) + delta = (x(k+1) - x(k)) / 2.0; + if (x(k+1) < x(k)) + error ("bar: x vector values must be in ascending order"); + endif + endif + k++; + endfor + else + error ("bar: arguments must be the same length"); + endif + else + error ("bar: arguments must be vectors"); + endif + else + error ("usage: [xb, yb] = bar (x, y)"); + endif + + if (nargout == 1) + plot (xb, yb); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/contour.m @@ -0,0 +1,62 @@ +function contour (z, n, x, y) + +# usage: contour (z, n, x, y) +# +# See also: plot, semilogx, semilogy, loglog, polar, mesh, contour, +# bar, stairs, gplot, gsplot, replot, xlabel, ylabel, title + + + if (nargin == 1) + n = 10; + endif + + if (nargin == 1 || nargin == 2) + if (is_matrix (z)) + set nosurface; + set contour; + set cntrparam bspline + command = sprintf ("set cntrparam levels %d", n); + eval (command); + set noparametric; + set view 0, 0, 1.9, 1 + gsplot z w l 1; + else + error ("mesh: argument must be a matrix"); + endif + elseif (nargin == 4) + if (is_vector (x) && is_vector (y) && is_matrix (z)) + xlen = length (x); + ylen = length (y); + if (xlen == rows (z) && ylen == columns (z)) + if (rows (x) == 1) + x = x'; + endif + len = 3 * ylen; + zz = zeros (xlen, ylen); + k = 1; + for i = 1:3:len + zz(:,i) = x; + zz(:,i+1) = y(k) * ones (xlen, 1); + zz(:,i+2) = z(:,k); + k++; + endfor + set nosurface + set contour + set cntrparam bspline + command = sprintf ("set cntrparam levels %d", n); + eval (command); + set parametric; + set view 0, 0, 1.9, 1 + gsplot zz w l 1; + else + disp ("mesh: rows (z) must be the same as length (x)"); + error (" and columns (z) must be the same as length (y)"); + endif + else + error ("mesh: x and y must be vectors and z must be a matrix"); + endif + else + error ("usage: mesh (z, levels, x, y)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/grid.m @@ -0,0 +1,30 @@ +function grid (x) + +# usage: grid ("on" | "off") +# +# Turn grid lines on or off for plotting. +# +# If the argument is omitted, "on" is assumed. +# +# See also: plot, semilogx, semilogy, loglog, polar, mesh, contour, +# bar, stairs, gplot, gsplot, replot, xlabel, ylabel, title + + if (nargin == 0) + set grid; + elseif (nargin == 1) + if (isstr (x)) + if (strcmp ("off", x)) + set nogrid; + elseif (strcmp ("on", x)) + set grid; + else + error ("usage: grid (\"on\" | \"off\")"); + endif + else + error ("error: grid: argument must be a string"); + endif + else + error ("usage: grid (\"on\" | \"off\")"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/loglog.m @@ -0,0 +1,28 @@ +function loglog (x1, x2) + +# usage: loglog (x, y) +# +# Make a 2D plot of y versus x using log scales for both axes. +# +# See the help message for the plot command for a description of how +# the arguments are interpreted. +# +# See also: plot, semilogx, semilogy, polar, mesh, contour, bar, stairs, +# gplot, gsplot, replot, xlabel, ylabel, title + + + set logscale x; + set logscale y; + set nopolar; + + if (nargin == 1) + plot_int (x1); + elseif (nargin == 2) + plot_int (x1, x2); + else + usage = sprintf ("usage: loglog (x)\n"); + usage = sprintf ("%s loglog (x, y)", usage); + error (usage); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/mesh.m @@ -0,0 +1,56 @@ +function mesh (x, y, z) + +# usage: mesh (x, y, z) +# +# See also: plot, semilogx, semilogy, loglog, polar, meshdom, contour, +# bar, stairs, gplot, gsplot, replot, xlabel, ylabel, title + + if (nargin == 1) + z = x; + if (is_matrix (z)) + set hidden3d; + set data style lines; + set surface; + set nocontour; + set noparametric; + set view 60, 30, 1, 1 + gsplot (z); + else + error ("mesh: argument must be a matrix"); + endif + elseif (nargin == 3) + if (is_vector (x) && is_vector (y) && is_matrix (z)) + xlen = length (x); + ylen = length (y); + if (xlen == rows (z) && ylen == columns (z)) + if (rows (x) == 1) + x = x'; + endif + len = 3 * ylen; + zz = zeros (xlen, ylen); + k = 1; + for i = 1:3:len + zz(:,i) = x; + zz(:,i+1) = y(k) * ones (xlen, 1); + zz(:,i+2) = z(:,k); + k++; + endfor + set hidden3d; + set data style lines; + set surface; + set nocontour; + set parametric; + set view 60, 30, 1, 1 + gsplot (zz); + else + disp ("mesh: rows (z) must be the same as length (x)"); + error (" and columns (z) must be the same as length (y)"); + endif + else + error ("mesh: x and y must be vectors and z must be a matrix"); + endif + else + error ("usage: mesh (z)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/meshdom.m @@ -0,0 +1,39 @@ +function [xx, yy] = meshdom (x, y) + +# usage: [xx, yy] = meshdom (x, y) +# +# Given vectors of x and y coordinates, return two matrices +# corresponding to the x and y coordinates of the mesh. +# +# See the file sombrero.m for an example of using mesh and meshdom. +# +# See also: plot, semilogx, semilogy, loglog, polar, mesh, contour, +# bar, stairs, gplot, gsplot, replot, xlabel, ylabel, title + + if (nargin == 2) + if (is_vector (x) && is_vector (y)) + xlen = length (x); + ylen = length (y); + xx = zeros (ylen, xlen); + yy = zeros (ylen, xlen); + y = y (ylen:-1:1); + if (columns (x) == 1) + x = x'; + endif + if (rows (y) == 1) + y = y'; + endif + for i = 1:ylen + xx(i,:) = x; + endfor + for i = 1:xlen + yy(:,i) = y; + endfor + else + error ("meshdom: arguments must be vectors"); + endif + else + error ("usage: [xx, yy] = meshdom (x, y)"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/plot.m @@ -0,0 +1,43 @@ +function plot (x1, x2) + +# usage: plot (x, y) +# +# If the first argument is a vector and the second is a matrix, the +# the vector is plotted versus the columns (or rows) of the matrix. +# (using whichever combination matches, with columns tried first.) +# +# If the first argument is a matrix and the second is a vector, the +# the columns (or rows) of the matrix are plotted versus the vector. +# (using whichever combination matches, with columns tried first.) +# +# If both arguments are vectors, the elements of y are plotted versus +# the elements of x. +# +# If both arguments are matrices, the columns of y are plotted versus +# the columns of x. In this case, both matrices must have the same +# number of rows and columns and no attempt is made to transpose the +# arguments to make the number of rows match. +# +# If both arguments are scalars, a single point is plotted. +# +# If only one argument is given, it is taken as the set of y +# coordinates and the x coordinates are taken to be the indices of the +# elements, starting with 1. +# +# See also: semilogx, semilogy, loglog, polar, mesh, contour, +# bar, stairs, gplot, gsplot, replot, xlabel, ylabel, title + + set nologscale; + set nopolar; + + if (nargin == 1) + plot_int (x1); + elseif (nargin == 2) + plot_int (x1, x2); + else + usage = sprintf ("usage: plot (x)\n"); + usage = sprintf ("%s plot (x, y)", usage); + error (usage); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/polar.m @@ -0,0 +1,23 @@ +function polar (x1, x2) + +# usage: polar (theta, rho) +# +# Make a 2D plot given polar the coordinates theta and rho. +# +# See also: plot, semilogx, semilogy, loglog, mesh, contour, bar, +# stairs, gplot, gsplot, replot, xlabel, ylabel, title + + set nologscale; + set nopolar; + + if (nargin == 1) + polar_int (x1); + elseif (nargin == 2) + polar_int (x1, x2); + else + usage = sprintf ("usage: polar (x)\n"); + usage = sprintf ("%s polar (x, y)", usage); + error (usage); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/semilogx.m @@ -0,0 +1,27 @@ +function semilogx (x1, x2) + +# usage: semilogx (x, y) +# +# Make a 2D plot of y versus x using a log scale for the x axis. +# +# See the help message for the plot command for a description of how +# the arguments are interpreted. +# +# See also: plot, semilogy, loglog, polar, mesh, contour, bar, stairs, +# gplot, gsplot, replot, xlabel, ylabel, title + + set logscale x; + set nologscale y; + set nopolar; + + if (nargin == 1) + plot_int (x1); + elseif (nargin == 2) + plot_int (x1, x2); + else + usage = sprintf ("usage: semilogx (x)\n"); + usage = sprintf ("%s semilogx (x, y)", usage); + error (usage); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/semilogy.m @@ -0,0 +1,27 @@ +function semilogy (x1, x2) + +# usage: semilogy (x, y) +# +# Make a 2D plot of y versus x using a log scale for the y axis. +# +# See the help message for the plot command for a description of how +# the arguments are interpreted. +# +# See also: plot, semilogx, loglog, polar, mesh, contour, bar, stairs, +# gplot, gsplot, replot, xlabel, ylabel, title + + set nologscale x; + set logscale y; + set nopolar; + + if (nargin == 1) + plot_int (x1); + elseif (nargin == 2) + plot_int (x1, x2); + else + usage = sprintf ("usage: semilogy (x)\n"); + usage = sprintf ("%s semilogy (x, y)", usage); + error (usage); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/sombrero.m @@ -0,0 +1,21 @@ +function sombrero (n) + +# usage: sombrero (n) +# +# Draw a `sombrero' in three dimensions using n grid lines. The +# function plotted is +# +# z = sin (x^2 + y^2) / (x^2 + y^2); + + if (nargin != 1) + error ("usage: sombrero (n)"); + endif + + x = y = linspace (-8, 8, n)'; + [xx, yy] = meshdom (x, y); + r = sqrt (xx .^ 2 + yy .^ 2) + eps; + z = sin (r) ./ r; + + mesh (x, y, z); + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/stairs.m @@ -0,0 +1,78 @@ +function [xs, ys] = stairs (x, y) + +# usage: [xs, ys] = bar (x, y) +# +# Given two vectors of x-y data, bar produces a `stairstep' plot. +# +# If only one argument is given, it is taken as a vector of y-values +# and the x coordiates are taken to be the indices of the elements. +# +# If two output arguments are specified, the data are generated but +# not plotted. For example, +# +# bar (x, y); +# +# and +# +# [xs, ys] = stairs (x, y); +# plot (xs, ys); +# +# are equivalent. +# +# See also: plot, semilogx, semilogy, loglog, polar, mesh, contour, +# bar, gplot, gsplot, replot, xlabel, ylabel, title + + + if (nargin == 1) + if (is_vector (x)) + len = 2 * length (x); + xs = ys = zeros (len, 1); + k = 0; + for i = 1:2:len + xs(i) = k++; + ys(i) = x(k); + ys(i+1) = x(k); + xs(i+1) = k; + endfor + else + error ("stairs: argument must be a vector"); + endif + elseif (nargin == 2) + if (is_vector (x) && is_vector (y)) + xlen = length (x); + ylen = length (y); + if (xlen == ylen) + len = 2 * xlen; + xs = ys = zeros (len, 1); + k = 1; + len_m2 = len - 2; + for i = 1:2:len_m2 + xs(i) = x(k); + ys(i) = y(k); + ys(i+1) = y(k); + k++; + xs(i+1) = x(k); + if (x(k) < x(k-1)) + error ("stairs: x vector values must be in ascending order"); + endif + endfor + xs(len-1) = x(xlen); + delta = x(xlen) - x(xlen-1); + xs(len) = x(xlen) + delta; + ys(len-1) = y(ylen); + ys(len) = y(ylen); + else + error ("stairs: arguments must be the same length"); + endif + else + error ("stairs: arguments must be vectors"); + endif + else + error ("usage: [xs, ys] = stairs (x, y)"); + endif + + if (nargout == 1) + plot (xs, ys); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/title.m @@ -0,0 +1,22 @@ +function title (text) + +# usage: title (text) +# +# Defines a title for a plot. The title will appear the next time a +# plot is displayed. +# +# See also: plot, semilogx, semilogy, loglog, polar, mesh, contour, +# bar, stairs, gplot, gsplot, replot, xlabel, ylabel + + if (nargin != 1) + error ("usage: title (text)"); + endif + + if (isstr (text)) + command = sprintf ("set title \"%s\"", text); + eval (command); + else + error ("error: title: text must be a string"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/xlabel.m @@ -0,0 +1,22 @@ +function xlabel (text) + +# usage: xlabel (text) +# +# Defines a label for the x-axis of a plot. The label will appear the +# next time a plot is displayed. +# +# See also: plot, semilogx, semilogy, loglog, polar, mesh, contour, +# bar, stairs, gplot, gsplot, replot, ylabel, title + + if (nargin != 1) + error ("usage: xlabel (text)"); + endif + + if (isstr (text)) + command = sprintf ("set xlabel \"%s\"", text); + eval (command); + else + error ("error: xlabel: text must be a string"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/plot/ylabel.m @@ -0,0 +1,22 @@ +function ylabel (text) + +# usage: ylabel (text) +# +# Defines a label for the y-axis of a plot. The label will appear the +# next time a plot is displayed. +# +# See also: plot, semilogx, semilogy, loglog, polar, mesh, contour, +# bar, stairs, gplot, gsplot, replot, xlabel, title + + if (nargin != 1) + error ("usage: ylabel (text)"); + endif + + if (isstr (text)) + command = sprintf ("set ylabel \"%s\"", text); + eval (command); + else + error ("error: ylabel: text must be a string"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/special-matrix/hadamard.m @@ -0,0 +1,20 @@ +function retval = hadamard (k) + +# usage: hadamard (k) +# +# Return the Hadamard matrix of order n = 2^k. +# +# See also: hankel, vander, hilb, invhilb, toeplitz + + if (nargin != 1) + error ("usage: hadamard (n)"); + endif + + if (k < 1) + retval = 1; + else + tmp = hadamard (k-1); + retval = [tmp, tmp; tmp, -tmp]; + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/special-matrix/hankel.m @@ -0,0 +1,59 @@ +function retval = hankel (c, r) + +# usage: hankel (c, r) +# +# Return the Hankel matrix constructed given the first column +# c, and (optionally) the last row r. +# +# If the second argument is omitted, the last row is taken to be the +# same as the first column. If the last element of c is not the same +# as the first element of r, the last element of c is used. +# +# See also: vander, hadamard, hilb, invhilb, toeplitz + + if (nargin == 1) + r = c; + elseif (nargin != 2) + error ("usage: hankel (c, r)"); + endif + + [c_nr, c_nc] = size (c); + [r_nr, r_nc] = size (r); + + if ((c_nr != 1 && c_nc != 1) || (r_nr != 1 && r_nc != 1)) + error ("hankel: expecting vector arguments") + endif + + if (c_nc != 1) + c = c'; + endif + + if (r_nr != 1) + r = r'; + endif + + if (r (1) != c (1)) + disp ("Column wins anti-diagonal conflict"); + endif + +# This should probably be done with the colon operator... + + nc = length (r); + nr = length (c); + + retval = zeros (nr, nc); + + for i = 1:min (nr, nc) + retval (1:nr-i+1, i) = c (i:nr); + endfor + + tmp = 1; + if (nc <= nr) + tmp = nr - nc + 2; + endif + + for i = nr:-1:tmp + retval (i, 2+nr-i:nc) = r (2:nc-nr+i); + endfor + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/special-matrix/hilb.m @@ -0,0 +1,29 @@ +function retval = hilb (n) + +# usage: hilb (n) +# +# Return the Hilbert matrix of order n. The i, j element of a Hilbert +# matrix is defined as +# +# H (i, j) = 1 / (i + j - 1); +# +# See also: hankel, vander, hadamard, invhilb, toeplitz + + + if (nargin != 1) + error ("usage: hilb (n)"); + endif + + nmax = length (n); + if (nmax == 1) + retval = zeros (n); + for j = 1:n + for i = 1:n + retval (i, j) = 1 / (i + j - 1); + endfor + endfor + else + error ("hilb: expecting scalar argument, found something else"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/special-matrix/invhilb.m @@ -0,0 +1,46 @@ +function retval = invhilb (n) + +# usage: invhilb (n) +# +# Return the inverse of a Hilbert matrix of order n. This is slow but +# exact. Compare with inv (hilb (n)). +# +# See also: hankel, vander, hadamard, hilb, toeplitz + + if (nargin != 1) + error ("usage: invhilb (n)"); + endif + + nmax = length (n); + if (nmax == 1) + retval = zeros (n); + for l = 1:n + for k = l:n + tmp = 1; + for i = 1:n + tmp = tmp * (i + k - 1); + endfor + for i = 1:n + if (i != k) + tmp = tmp * (l + i - 1); + endif + endfor + for i = 1:n + if (i != l) + tmp = tmp / (i - l); + endif + endfor + for i = 1:n + if (i != k) + tmp = tmp / (i - k); + endif + endfor + retval (k, l) = tmp; + retval (l, k) = tmp; + endfor + endfor + else + error ("hilb: expecting scalar argument, found something else"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/special-matrix/toeplitz.m @@ -0,0 +1,54 @@ +function retval = toeplitz (c, r) + +# usage: toeplitz (c, r) +# +# Return the Toeplitz matrix constructed given the first column +# c, and (optionally) the first row r. +# +# If the second argument is omitted, the first row is taken to be the +# same as the first column. If the first element of c is not the same +# as the first element of r, the first element of c is used. +# +# See also: hankel, vander, hadamard, hilb, invhib + + if (nargin == 1) + r = c; + elseif (nargin != 2) + error ("usage: toeplitz (c, r)"); + endif + + [c_nr, c_nc] = size (c); + [r_nr, r_nc] = size (r); + + if ((c_nr != 1 && c_nc != 1) || (r_nr != 1 && r_nc != 1)) + error ("toeplitz: expecting vector arguments") + endif + + if (c_nc != 1) + c = c'; + endif + + if (r_nr != 1) + r = r'; + endif + + if (r (1) != c (1)) + disp ("Column wins diagonal conflict"); + endif + +# This should probably be done with the colon operator... + + nc = length (r); + nr = length (c); + + retval = zeros (nr, nc); + + for i = 1:min (nc, nr) + retval (i:nr, i) = c (1:nr-i+1); + endfor + + for i = 1:min (nr, nc-1) + retval (i, i+1:nc) = r (2:nc-i+1); + endfor + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/special-matrix/vander.m @@ -0,0 +1,32 @@ +function retval = vander (c) + +# usage: vander (c) +# +# Return the Vandermonde matrix whose next to last column is c. +# +# See also: hankel, hadamard, hilb, invhilb, toeplitz + + if (nargin != 1) + error ("usage: vander (c)"); + endif + + nr = rows (c); + nc = columns (c); + if (nr == 1 && nc == 1) + retval = 1; + elseif (nr == 1 || nc == 1) + n = length (c); + if (n > 0) + retval = zeros (n, n); + for i = 1:n + tmp = c(i); + for j = 1:n + retval (i, j) = tmp ^ (n - j); + endfor + endfor + endif + else + error ("vander: argument must be a vector"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/statistics/mean.m @@ -0,0 +1,25 @@ +function retval = mean (a) + +# usage: mean (a) +# +# For vector arguments, return the mean the values. +# +# For matrix arguments, return a row vector containing the mean for +# each column. +# +# See also: median, std + + if (nargin != 1) + error ("usage: mean (a)"); + endif + + [nr, nc] = size (a); + if (nr == 1 || nc == 1) + retval = sum (a) / length (a); + elseif (nr > 0 && nc > 0) + retval = sum (a) / nr; + else + error ("mean: invalid matrix argument"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/statistics/median.m @@ -0,0 +1,38 @@ +function retval = median (a) + +# usage: median (a) +# +# For vector arguments, return the median of the values. +# +# For matrix arguments, return a row vector containing the median for +# each column. +# +# See also: std, mean + + if (nargin != 1) + error ("usage: medain (a)"); + endif + + [nr, nc] = size (a); + s = sort (a); + if (nr == 1 && nc > 0) + if (rem (nc, 2) == 0) + i = nc/2; + retval = (s (i) + s (i+1)) / 2; + else + i = ceil (nc/2); + retval = s (i); + endif + elseif (nr > 0 && nc > 0) + if (rem (nr, 2) == 0) + i = nr/2; + retval = (s (i,:) + s (i+1,:)) / 2; + else + i = ceil (nr/2); + retval = s (i,:); + endif + else + error ("median: invalid matrix argument"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/statistics/std.m @@ -0,0 +1,30 @@ +function retval = std (a) + +# usage: std (a) +# +# For vector arguments, std returns the standard deviation of the +# values. For matrix arguments, std returns a row vector containing +# the standard deviation for each column. +# +# See also: mean, median + + if (nargin != 1) + error ("usage: std (a)"); + endif + + nr = rows (a); + nc = columns (a); + if (nc == 1 && nr == 1) + retval = 0; + elseif (nc == 1 || nr == 1) + tmp = sum (a); + n = length (a); + retval = sqrt ((n * sumsq (a) - tmp .* tmp) / (n * (n - 1))); + elseif (nr > 1 && nc > 0) + tmp = sum (a); + retval = sqrt ((nr * sumsq (a) - tmp .* tmp) / (nr * (nr - 1))); + else + error ("mean: invalid matrix argument"); + endif + +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/strings/strcmp.m @@ -0,0 +1,23 @@ +function status = strcmp (s1, s2) + +# usage: strcmp (s1, s2) +# +# Compare two strings. +# +# WARNING: Unlike the C function of the same name, this function +# returns 1 for equal and zero for not equal. Why? To be compatible +# with Matlab, of course. + + if (nargin != 2) + error ("usage: strcmp (s, t)"); + endif + + status = 0; + if (isstr (s1) && isstr(s2) && length (s1) == length (s2)) + tmp = implicit_str_to_num_ok; + implicit_str_to_num_ok = "true"; + status = all (s1 == s2); + implicit_str_to_num_ok = tmp; + endif + +endfunction