text:
function test_inverse ( )

%*****************************************************************************80
%
%% TEST_INVERSE tests the inverse computations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_INVERSE\n' );
  fprintf ( 1, '  A = a test matrix of order N;\n' );
  fprintf ( 1, '  B = inverse as computed by a routine.\n' );
  fprintf ( 1, '  C = inverse as computed by Matlab''s INV function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ||A||    = Frobenius norm of A.\n' );
  fprintf ( 1, '  ||C||    = Frobenius norm of C.\n' );
  fprintf ( 1, '  ||I-AC|| = Frobenius norm of I-A*C.\n' );
  fprintf ( 1, '  ||I-AB|| = Frobenius norm of I-A*B.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Title                    N      ' );
  fprintf ( 1, '||A||      ||C||      ||I-AC||        ||I-AB||\n' );
  fprintf ( 1, '\n' );
%
%  AEGERTER
%
  title = 'AEGERTER';
  n = 5;
  a = aegerter ( n );
  b = aegerter_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  BAB
%
  title = 'BAB';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  [ beta, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = bab ( n, alpha, beta );
  b = bab_inverse ( n, alpha, beta );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  BAUER
%
  title = 'BAUER';
  n = 6;
  a = bauer ( );
  b = bauer_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  BERNSTEIN
%
  title = 'BERNSTEIN';
  n = 5;
  a = bernstein ( n );
  b = bernstein_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  BIS
%
  title = 'BIS';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  [ beta, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = bis ( alpha, beta, n, n );
  b = bis_inverse ( alpha, beta, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  BIW
%
  title = 'BIW';
  n = 5;
  a = biw ( n );
  b = biw_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  BODEWIG
%
  title = 'BODEWIG';
  n = 4;
  a = bodewig ( );
  b = bodewig_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  BOOTHROYD
%
  title = 'BOOTHROYD';
  n = 5;
  a = boothroyd ( n );
  b = boothroyd_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  BORDERBAND
%
  title = 'BORDERBAND';
  n = 5;
  a = borderband ( n );
  b = borderband_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CARRY
%
  title = 'CARRY';
  n = 5;
  i4_lo = 2;
  i4_hi = 20;
  seed = 123456789;
  [ k, seed ] = i4_uniform_ab ( i4_lo, i4_hi, seed );
  a = carry ( n, k );
  b = carry_inverse ( n, k );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CAUCHY
%
  title = 'CAUCHY';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  [ y, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  a = cauchy ( n, x, y );
  b = cauchy_inverse ( n, x, y );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CHEBY_T
%
  title = 'CHEBY_T';
  n = 5;
  a = cheby_t ( n );
  b = cheby_t_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CHEBY_U
%
  title = 'CHEBY_U';
  n = 5;
  a = cheby_u ( n );
  b = cheby_u_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CHEBY_VAN2
%
  title = 'CHEBY_VAN2';
  n = 5;
  a = cheby_van2 ( n );
  b = cheby_van2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CHEBY_VAN3
%
  title = 'CHEBY_VAN3';
  n = 5;
  a = cheby_van3 ( n );
  b = cheby_van3_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CHOW
%
  title = 'CHOW';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  [ beta, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = chow ( alpha, beta, n, n );
  b = chow_inverse ( alpha, beta, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CIRCULANT
%
  title = 'CIRCULANT';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  a = circulant ( n, n, x );
  b = circulant_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CIRCULANT2
%
  title = 'CIRCULANT2';
  n = 5;
  a = circulant2 ( n );
  b = circulant2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CLEMENT1
%  N must be even.
%
  title = 'CLEMENT1';
  n = 6;
  a = clement1 ( n );
  b = clement1_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CLEMENT2
%  N must be even.
%
  title = 'CLEMENT2';
  n = 6;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n - 1, r8_lo, r8_hi, seed );
  [ y, seed ] = r8vec_uniform_ab ( n - 1, r8_lo, r8_hi, seed );
  a = clement2 ( n, x, y );
  b = clement2_inverse ( n, x, y );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  COMBIN
%
  title = 'COMBIN';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  [ beta, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = combin ( alpha, beta, n );
  b = combin_inverse ( alpha, beta, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  COMPANION
%
  title = 'COMPANION';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  a = companion ( n, x );
  b = companion_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  COMPLEX_I
%
  title = 'COMPLEX_I';
  n = 2;
  a = complex_i ( );
  b = complex_i_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CONEX1
%
  title = 'CONEX1';
  n = 4;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = conex1 ( alpha );
  b = conex1_inverse ( alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CONEX2
%
  title = 'CONEX2';
  n = 3;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = conex2 ( alpha );
  b = conex2_inverse ( alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CONEX3
%
  title = 'CONEX3';
  n = 5;
  a = conex3 ( n );
  b = conex3_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  CONFERENCE
%
  title = 'CONFERENCE';
  n = 6;
  a = conference ( n );
  b = conference_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DAUB2
%
  title = 'DAUB2';
  n = 4;
  a = daub2 ( n );
  b = daub2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DAUB4
%
  title = 'DAUB4';
  n = 8;
  a = daub4 ( n );
  b = daub4_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DAUB6
%
  title = 'DAUB6';
  n = 12;
  a = daub6 ( n );
  b = daub6_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DAUB8
%
  title = 'DAUB8';
  n = 16;
  a = daub8 ( n );
  b = daub8_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DAUB10
%
  title = 'DAUB10';
  n = 20;
  a = daub10 ( n );
  b = daub10_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DAUB12
%
  title = 'DAUB12';
  n = 24;
  a = daub12 ( n );
  b = daub12_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DIAGONAL
%
  title = 'DIAGONAL';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  a = diagonal ( n, n, x );
  b = diagonal_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DIF1
%  N must be even.
%
  title = 'DIF1';
  n = 6;
  a = dif1 ( n, n );
  b = dif1_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DIF2
%
  title = 'DIF2';
  n = 5;
  a = dif2 ( n, n );
  b = dif2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DORR
%
  title = 'DORR';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = dorr ( alpha, n );
  b = dorr_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  DOWNSHIFT
%
  title = 'DOWNSHIFT';
  n = 5;
  a = downshift ( n );
  b = downshift_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  EULERIAN
%
  title = 'EULERIAN';
  n = 5;
  a = eulerian ( n, n );
  b = eulerian_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  EXCHANGE
%
  title = 'EXCHANGE';
  n = 5;
  a = exchange ( n, n );
  b = exchange_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  FIBONACCI2
%
  title = 'FIBONACCI2';
  n = 5;
  a = fibonacci2 ( n );
  b = fibonacci2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  FIBONACCI3
%
  title = 'FIBONACCI3';
  n = 5;
  a = fibonacci3 ( n );
  b = fibonacci3_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  FIEDLER.
%  The FIEDLER_INVERSE routine assumes the X vector is sorted.
%
  title = 'FIEDLER';
  n = 7;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  x = sort ( x );
  a = fiedler ( n, n, x );
  b = fiedler_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  FORSYTHE
%
  title = 'FORSYTHE';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  [ beta, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = forsythe ( alpha, beta, n );
  b = forsythe_inverse ( alpha, beta, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  FOURIER_COSINE
%
  title = 'FOURIER_COSINE';
  n = 5;
  a = fourier_cosine ( n );
  b = fourier_cosine_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  FOURIER_SINE
%
  title = 'FOURIER_SINE';
  n = 5;
  a = fourier_sine ( n );
  b = fourier_sine_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  FRANK
%
  title = 'FRANK';
  n = 5;
  a = frank ( n );
  b = frank_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  GFPP
%
  title = 'GFPP';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = gfpp ( n, alpha );
  b = gfpp_inverse ( n, alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  GIVENS
%
  title = 'GIVENS';
  n = 5;
  a = givens ( n, n );
  b = givens_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  GK316
%
  title = 'GK316';
  n = 5;
  a = gk316 ( n );
  b = gk316_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  GK323
%
  title = 'GK323';
  n = 5;
  a = gk323 ( n, n );
  b = gk323_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  GK324
%
  title = 'GK324';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n - 1, r8_lo, r8_hi, seed );
  a = gk324 ( n, n, x );
  b = gk324_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  HANKEL_N
%
  title = 'HANKEL_N';
  n = 6;
  a = hankel_n ( n );
  b = hankel_n_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  HANOWA
%
  title = 'HANOWA';
  n = 6;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = hanowa ( alpha, n );
  b = hanowa_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  HARMAN
%
  title = 'HARMAN';
  n = 8;
  a = harman ( );
  b = harman_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  HARTLEY
%
  title = 'HARTLEY';
  n = 5;
  a = hartley ( n );
  b = hartley_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  HELMERT
%
  title = 'HELMERT';
  n = 5;
  a = helmert ( n );
  b = helmert_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  HELMERT2
%
  title = 'HELMERT2';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  a = helmert2 ( n, x );
  b = helmert2_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  HERMITE
%
  title = 'HERMITE';
  n = 5;
  a = hermite ( n );
  b = hermite_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  HERNDON
%
  title = 'HERNDON';
  n = 5;
  a = herndon ( n );
  b = herndon_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  HILBERT
%
  title = 'HILBERT';
  n = 5;
  a = hilbert ( n, n );
  b = hilbert_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  HOUSEHOLDER
%
  title = 'HOUSEHOLDER';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  a = householder ( n, x );
  b = householder_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  IDENTITY
%
  title = 'IDENTITY';
  n = 5;
  a = identity ( n, n );
  b = identity_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  ILL3
%
  title = 'ILL3';
  n = 3;
  a = ill3 ( );
  b = ill3_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  INTEGRATION
%
  title = 'INTEGRATION';
  n = 6;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = integration ( alpha, n );
  b = integration_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  INVOL
%
  title = 'INVOL';
  n = 5;
  a = invol ( n );
  b = invol_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  JACOBI
%  N must be even.
%
  title = 'JACOBI';
  n = 6
  a = jacobi ( n, n );
  b = jacobi_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  JORDAN
%
  title = 'JORDAN';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = jordan ( n, n, alpha );
  b = jordan_inverse ( n, alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  KAHAN
%
  title = 'KAHAN';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = kahan ( alpha, n, n );
  b = kahan_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  KERSHAW
%
  title = 'KERSHAW';
  n = 4;
  a = kershaw ( );
  b = kershaw_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  KERSHAWTRI
%
  title = 'KERSHAWTRI';
  n = 5;
  x_n = floor ( ( n + 1 ) / 2 );
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( x_n, r8_lo, r8_hi, seed );
  a = kershawtri ( n, x );
  b = kershawtri_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  KMS
%
  title = 'KMS';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = kms ( alpha, n, n );
  b = kms_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  LAGUERRE
%
  title = 'LAGUERRE';
  n = 5;
  a = laguerre ( n );
  b = laguerre_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  LEGENDRE
%
  title = 'LEGENDRE';
  n = 5;
  a = legendre ( n );
  b = legendre_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  LEHMER
%
  title = 'LEHMER';
  n = 5;
  a = lehmer ( n, n );
  b = lehmer_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  LESP
%
  title = 'LESP';
  n = 5;
  a = lesp ( n, n );
  b = lesp_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  LIETZKE
%
  title = 'LIETZKE';
  n = 5;
  a = lietzke ( n );
  b = lietzke_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  LINE_ADJ
%  N must be even.
%
  title = 'LINE_ADJ';
  n = 6;
  a = line_adj ( n );
  b = line_adj_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  LOTKIN
%
  title = 'LOTKIN';
  n = 5;
  a = lotkin ( n, n );
  b = lotkin_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  MAXIJ
%
  title = 'MAXIJ';
  n = 5;
  a = maxij ( n, n );
  b = maxij_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  MILNES
%
  title = 'MILNES';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  a = milnes ( n, n, x );
  b = milnes_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  MINIJ
%
  title = 'MINIJ';
  n = 5;
  a = minij ( n, n );
  b = minij_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  MOLER1
%
  title = 'MOLER1';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = moler1 ( alpha, n, n );
  b = moler1_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  MOLER3
%
  title = 'MOLER3';
  n = 5;
  a = moler3 ( n, n );
  b = moler3_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  ORTEGA
%
  title = 'ORTEGA';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ v1, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  [ v2, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  [ v3, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  a = ortega ( n, v1, v2, v3 );
  b = ortega_inverse ( n, v1, v2, v3 );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  ORTH_SYMM
%
  title = 'ORTH_SYMM';
  n = 5;
  a = orth_symm ( n );
  b = orth_symm_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  OTO
%
  title = 'OTO';
  n = 5;
  a = oto ( n, n );
  b = oto_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  PARTER
%
  title = 'PARTER';
  n = 5;
  a = parter ( n, n );
  b = parter_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  PASCAL1
%
  title = 'PASCAL1';
  n = 5;
  a = pascal1 ( n );
  b = pascal1_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  PASCAL2
%
  title = 'PASCAL2';
  n = 5;
  a = pascal2 ( n );
  b = pascal2_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  PASCAL3
%
  title = 'PASCAL3';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = pascal3 ( n, alpha );
  b = pascal3_inverse ( n, alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  PDS_RANDOM
%
  title = 'PDS_RANDOM';
  n = 5;
  key = 123456789;
  a = pds_random ( n, key );
  b = pds_random_inverse ( n, key );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  PEI
%
  title = 'PEI';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = pei ( alpha, n );
  b = pei_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  PERMUTATION_RANDOM
%
  title = 'PERMUTATION_RANDOM';
  n = 5;
  key = 123456789;
  a = permutation_random ( n, key );
  b = permutation_random_inverse ( n, key );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  PLU
%
  title = 'PLU';
  n = 5;
  seed = 123456789;
  for i = 1 : n
    i4_lo = i;
    i4_hi = n;
    [ j, seed ] = i4_uniform_ab ( i4_lo, i4_hi, seed );
    pivot(i) = j;
  end
  a = plu ( n, pivot );
  b = plu_inverse ( n, pivot );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  RIS
%
  title = 'RIS';
  n = 5;
  a = ris ( n );
  b = ris_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  RODMAN
%
  title = 'RODMAN';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = rodman ( n, n, alpha );
  b = rodman_inverse ( n, alpha );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  RUTIS1
%
  title = 'RUTIS1';
  n = 4;
  a = rutis1 ( );
  b = rutis1_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  RUTIS2
%
  title = 'RUTIS2';
  n = 4;
  a = rutis2 ( );
  b = rutis2_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  RUTIS3
%
  title = 'RUTIS3';
  n = 4;
  a = rutis3 ( );
  b = rutis3_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  RUTIS4
%
  title = 'RUTIS4';
  n = 4;
  a = rutis4 ( n );
  b = rutis4_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  RUTIS5
%
  title = 'RUTIS5';
  n = 4;
  a = rutis5 ( );
  b = rutis5_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  SCHUR_BLOCK
%
  title = 'SCHUR_BLOCK';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  x_n = floor ( ( n + 1 ) / 2 );
  [ x, seed ] = r8vec_uniform_ab ( x_n, r8_lo, r8_hi, seed );
  y_n = floor ( n / 2 );
  [ y, seed ] = r8vec_uniform_ab ( y_n, r8_lo, r8_hi, seed );
  a = schur_block ( n, x, y );
  b = schur_block_inverse ( n, x, y );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  SPLINE
%
  title = 'SPLINE';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n - 1, r8_lo, r8_hi, seed );
  a = spline ( n, x );
  b = spline_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  STIRLING
%
  title = 'STIRLING';
  n = 5;
  a = stirling ( n, n );
  b = stirling_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  SUMMATION
%
  title = 'SUMMATION';
  n = 5;
  a = summation ( n, n );
  b = summation_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  SWEET1
%
  title = 'SWEET1';
  n = 6;
  a = sweet1 ( );
  b = sweet1_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  SWEET2
%
  title = 'SWEET2';
  n = 6;
  a = sweet2 ( );
  b = sweet2_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  SWEET3
%
  title = 'SWEET3';
  n = 6;
  a = sweet3 ( );
  b = sweet3_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  SWEET4
%
  title = 'SWEET4';
  n = 13;
  a = sweet4 ( );
  b = sweet4_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  SYLVESTER_KAC
%  N must be even.
%
  title = 'SYLVESTER_KAC';
  n = 6;
  a = sylvester_kac ( n );
  b = sylvester_kac_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  SYMM_RANDOM
%
  title = 'SYMM_RANDOM';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  key = 123456789;
  a = symm_random ( n, x, key );
  b = symm_random_inverse ( n, x, key );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  TRI_UPPER
%
  title = 'TRI_UPPER';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = tri_upper ( alpha, n );
  b = tri_upper_inverse ( alpha, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  TRIS
%
  title = 'TRIS';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  [ beta, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  [ gamma, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = tris ( n, n, alpha, beta, gamma );
  b = tris_inverse ( n, alpha, beta, gamma );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  TRIV
%
  title = 'TRIV';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n - 1, r8_lo, r8_hi, seed );
  [ y, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  [ z, seed ] = r8vec_uniform_ab ( n - 1, r8_lo, r8_hi, seed );
  a = triv ( n, x, y, z );
  b = triv_inverse ( n, x, y, z );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  TRIW
%
  title = 'TRIW';
  n = 5;
  i4_lo = 0;
  i4_hi = n - 1;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ k, seed ] = i4_uniform_ab ( i4_lo, i4_hi, seed );
  [ alpha, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
  a = triw ( alpha, k, n );
  b = triw_inverse ( alpha, k, n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  UNITARY_RANDOM
%
  title = 'UNITARY_RANDOM';
  n = 5;
  key = 123456789;
  a = unitary_random ( n, key );
  b = unitary_random_inverse ( n, key );
  c = inv ( a );
  error_ab = c8mat_is_inverse ( n, a, b );
  error_ac = c8mat_is_inverse ( n, a, c );
  norma_frobenius = c8mat_norm_fro ( n, n, a );
  normc_frobenius = c8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  UPSHIFT
%
  title = 'UPSHIFT';
  n = 5;
  a = upshift ( n );
  b = upshift_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  VAND1
%
  title = 'VAND1';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  a = vand1 ( n, x );
  b = vand1_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  VAND2
%
  title = 'VAND2';
  n = 5;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );
  a = vand2 ( n, x );
  b = vand2_inverse ( n, x );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  WILK03
%
  title = 'WILK03';
  n = 3;
  a = wilk03 ( );
  b = wilk03_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  WILK04
%
  title = 'WILK04';
  n = 4;
  a = wilk04 ( );
  b = wilk04_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  WILK05
%
  title = 'WILK05';
  n = 5;
  a = wilk05 ( );
  b = wilk05_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  WILK21
%
  title = 'WILK21';
  n = 21;
  a = wilk21 ( n );
  b = wilk21_inverse ( n );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );
%
%  WILSON
%
  title = 'WILSON';
  n = 4;
  a = wilson ( );
  b = wilson_inverse ( );
  c = inv ( a );
  error_ab = r8mat_is_inverse ( n, a, b );
  error_ac = r8mat_is_inverse ( n, a, c );
  norma_frobenius = r8mat_norm_fro ( n, n, a );
  normc_frobenius = r8mat_norm_fro ( n, n, c );
  fprintf ( 1, '  %-20s  %4d  %10.2g  %10.2g  %14g  %14g\n', ...
    title, n, norma_frobenius, normc_frobenius, error_ac, error_ab );

  return
end


text:
%% This function will exclude the left hand side guess from the right hand side library.
% Last Updated: 2019/04/22
% Coded By: K
function [RHS_Data,RHS_Struct]=ExcludeGuess(SINDy_Data,SINDy_Struct,LHS_Sym)
% First get the size of the library
[n,m]=size(SINDy_Data);
% Then run a for loop and generate a new library
j=1;
for i=1:m
    if SINDy_Struct{i}~=LHS_Sym
        RHS_Data(:,j)=SINDy_Data(:,i);
        RHS_Struct(1,j)=SINDy_Struct(1,i);
        j=j+1;
    else
        % For use of debuging
        SINDy_Struct{i};
    end
end

author: dynamicslab

repo: SINDy-PI


text:
function val=intPowSinPow(u,n,m)
%%INTPOWSINPOW Evaluate the integral of u^n*sin(u)^m du.  A definite
%           integral can be evaluated, or an indefinite integral (with a
%           particular additive constant).
%
%INPUTS: u A 2XN (for definite integral) or a 1XN (for indefinite
%          integrals) set of N points. For definite integrals, u(1,:) are
%          the real lower bounds and u(2,:) are the real upper bounds. For
%          indefinite integrals, the integral is evaluated at the points in
%          u. The values in u should be real.
%        n The positive integer exponent of u.
%        m The positive integer exponent of the sine term.
%
%OUTPUTS: val The 1XN set of values of the integral of u^n*sin(u)^m.
%
%This function implements formulas 4 and 5 of Section 2.631 of [1].
%
%REFERENCES:
%[1] I. S. Gradshteyn and I. M. Ryzhik, Tables of Integrals, Series, and
%    Products, Corrected and Enlarged Edition. New York: Academic Press,
%    1980, translated from Russian.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(u,1);

if(isempty(u))
   val=[];
   return;
end

if(numDim==1)%An indefinite integral
    val=indefIntPowSinPow(u,n,m);
else%A definite integral
    val=indefIntPowSinPow(u(2,:),n,m)-indefIntPowSinPow(u(1,:),n,m);
end

end

function val=indefIntPowSinPow(u,n,m)

if(mod(m,2)==0)%Use formula 4.
    m=m/2;

    val=0;
    for k=0:(m-1)
        a=2*(m-k);
        x=a*u;
        val=val+(-1)^k*binomial(2*m,k)*(1/a)^(n+1)*intPowCos(x,n);
    end

    val=(-1)^m/(2^(2*m-1))*val;
    val=val+binomial(2*m,m)*u.^(n+1)/(2^(2*m)*(n+1));
else%Use formula 5.
    m=(m-1)/2;

    val=0;
    for k=0:m
        a=2*(m-k)+1;
        x=a*u;
        val=val+(-1)^k*binomial(2*m+1,k)*(1/a)^(n+1)*intPowSin(x,n);
    end

    val=val*(-1)^m/(2^(2*m));
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [p] = gui_eeg_contour_steps(parent)

% gui_eeg_contour_steps - GUI for contour steps of topographic maps
%
% Usage: [p] = eeg_contour_steps([parent])
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:56 $

% Licence:  GNU GPL, no express or implied warranties
% History:  10/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% reset the gui_struct parent field to the handles of the input gui_struct

if exist('parent','var'),
    ContourSteps.parent.gui = parent;
    ParentUserdata = get(parent,'Userdata');
    ContourSteps.p = ParentUserdata.p;
else
    error('No parent gui parameter supplied');
end

% GUI General Parameters
GUIwidth  = 200;
GUIheight = 100;
GUI = figure('Name','Contour Steps','Tag','CONTOURSTEPS',...
             'NumberTitle','off','HandleVisibility','callback',...
             'MenuBar','none','Position',[1 1 GUIwidth GUIheight]);
movegui(GUI,'center');

ContourSteps.gui = GUI;

Font.FontName   = 'Helvetica';
Font.FontUnits  = 'Pixels';
Font.FontSize   = 12;
Font.FontWeight = 'normal';
Font.FontAngle  = 'normal';

% Contour Step Size

switch ContourSteps.p.contour.stepMethod
case 1,     step = [0 1];
otherwise,  step = [1 0];
end
% Step Size
G.BstepMethod0 = uicontrol('Parent',GUI,'Style','Radiobutton',...
    'Units','Normalized',Font,'Position',[.1 .6 .6 .2],...
    'String','Step Size','Value',step(1),...
    'Callback',strcat('ContourSteps = get(gcbf,''Userdata'');',...
                      'ContourSteps.p.contour.stepMethod = 0;',...
                      'set(ContourSteps.handles.BstepMethod0,''Value'',1);',...
                      'set(ContourSteps.handles.BstepMethod1,''Value'',0);',...
                      'set(ContourSteps.handles.EstepMethod0,''String'',num2str(ContourSteps.p.contour.stepSize));',...
                      'set(ContourSteps.gui,''Userdata'',ContourSteps); clear ContourSteps;'));
G.EstepMethod0 = uicontrol('Parent',GUI,'Style','edit',...
    'Units','Normalized',Font,'Position',[.7 .6 .2 .2],...
    'String',num2str(ContourSteps.p.contour.stepSize),...
    'Callback',strcat('ContourSteps = get(gcbf,''Userdata'');',...
                      'ContourSteps.p.contour.stepSize = str2num(get(ContourSteps.handles.EstepMethod0,''String''));',...
                      'set(ContourSteps.gui,''Userdata'',ContourSteps); clear ContourSteps;'));

% Number of steps
G.BstepMethod1 = uicontrol('Parent',GUI,'Style','Radiobutton',...
    'Units','Normalized',Font,'Position',[.1 .4 .6 .2],...
    'String','Number of Steps','Value',step(2),...
    'Callback',strcat('ContourSteps = get(gcbf,''Userdata'');',...
                      'ContourSteps.p.contour.stepMethod = 1;',...
                      'set(ContourSteps.handles.BstepMethod0,''Value'',0);',...
                      'set(ContourSteps.handles.BstepMethod1,''Value'',1);',...
                      'set(ContourSteps.handles.EstepMethod1,''String'',num2str(ContourSteps.p.contour.Nsteps));',...
                      'set(ContourSteps.gui,''Userdata'',ContourSteps); clear ContourSteps;'));
G.EstepMethod1 = uicontrol('Parent',GUI,'Style','edit',...
    'Units','Normalized',Font,'Position',[.7 .4 .2 .2],...
    'String',num2str(ContourSteps.p.contour.Nsteps),...
    'Callback',strcat('ContourSteps = get(gcbf,''Userdata'');',...
                      'ContourSteps.p.contour.Nsteps = str2num(get(ContourSteps.handles.EstepMethod1,''String''));',...
                      'set(ContourSteps.gui,''Userdata'',ContourSteps); clear ContourSteps;'));

Font.FontWeight = 'bold';

% OK: Return the parameters!
G.Bdone = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
    'Position',[.3 .15 .3 .2],...
    'String','OK','BusyAction','queue',...
    'TooltipString','Update the ascii parameters.',...
    'BackgroundColor',[0.0 0.5 0.0],...
    'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
    'Callback',strcat('ContourSteps = get(gcbf,''Userdata'');',...
                      'if isfield(ContourSteps,''parent''),',...
                          'if isfield(ContourSteps.parent,''gui''),',...
                            'parent = get(ContourSteps.parent.gui,''UserData'');',...
                            'if isfield(parent,''p'') & isfield(ContourSteps,''p''),',...
                                'parent.p = ContourSteps.p;',...
                                'set(ContourSteps.parent.gui,''UserData'',parent);',...
                            'end; ',...
                          'end; ',...
                      'end; ',...
                      'close gcbf;[p] = ContourSteps.p; clear ContourSteps parent;'));

% Cancel
G.Bquit = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
    'Position',[.6 .15 .3 .2],...
    'String','CANCEL','BusyAction','queue',...
    'BackgroundColor',[0.75 0.0 0.0],...
    'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
    'Callback','close gcbf;');

% Update the gui_struct handles for this gui
ContourSteps.handles = G;
set(ContourSteps.gui,'Userdata',ContourSteps);

p = ContourSteps.p;
return

author: PatternRecognition

repo: OpenBMI


text:
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
%   [all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
%   logisitc regression classifiers and returns each of these classifiers
%   in a matrix all_theta, where the i-th row of all_theta corresponds
%   to the classifier for label i

% Some useful variables
m = size(X, 1);
n = size(X, 2);

% You need to return the following variables correctly
all_theta = zeros(num_labels, n + 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
%               logistic regression classifiers with regularization
%               parameter lambda.
%
% Hint: theta(:) will return a column vector.
%
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell use
%       whether the ground truth is true/false for this class.
%
% Note: For this assignment, we recommend using fmincg to optimize the cost
%       function. It is okay to use a for-loop (for c = 1:num_labels) to
%       loop over the different classes.
%
%       fmincg works similarly to fminunc, but is more efficient when we
%       are dealing with large number of parameters.
%
% Example Code for fmincg:
%
%     % Set Initial theta
%     initial_theta = zeros(n + 1, 1);
%
%     % Set options for fminunc
%     options = optimset('GradObj', 'on', 'MaxIter', 50);
%
%     % Run fmincg to obtain the optimal theta
%     % This function will return theta and the cost
%     [theta] = ...
%         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
%                 initial_theta, options);
%

for c = 1:num_labels

    initial_theta = zeros(n + 1, 1);

    options = optimset('GradObj', 'on', 'MaxIter', 50);

    [theta] = ...
         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
                 initial_theta, options);

    all_theta(c,:) = theta';

end

% =========================================================================

end

author: rmarquis

repo: coursera-machinelearning


text:
function laguerre_test_int_test05 ( )

%*****************************************************************************80
%
%% LAGUERRE_TEST_INT_TEST05 tests P00_RAT_TRANSFORM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 December 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGUERRE_TEST_INT_TEST05\n' );
  fprintf ( 1, '  P00_RAT_TRANSFORM. applies a rational tranform\n' );
  fprintf ( 1, '  to estimate an integral on [ALPHA,+oo)\n' );
  fprintf ( 1, '  as a transformed integral on (0,1/(1+ALPHA)]\n' );
  fprintf ( 1, '  and applying a Gauss-Legendre rule.\n' );

  problem_num = p00_problem_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '                           Exact\n' );
  fprintf ( 1, '   Problem     Order       Estimate    Error\n' );

  for problem = 1 : problem_num

    exact = p00_exact ( problem );

    order = 1;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %8d            %14.6f\n', problem, exact );

    for order_log = 0 : 6

      estimate = p00_rat_transform ( problem, order );

      err = abs ( exact - estimate );

      fprintf ( 1, '            %8d  %14.6f  %14.6e\n', order, estimate, err );

      order = order * 2;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [accu,prediction] = kernel_knn_classification(test_kernel,train_label,K,test_label)
% this function performs classification with given kernels and labels
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input parameters:
% test_kernel: the test kernel
% train_label: the training label
% K: Ks of knn classifier
% test_label: the test label
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Output parameters:
% accu: the classification accuracy
% prediction: the predicted label for the test data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Jianjia Zhang, jz163@uowmail.edu.au Dec, 2014, all rights reserved
% For implementation details, please refer to:
% "Learning Discriminative Stein Kernel for SPD Matrices and Its Applications."
% arXiv preprint arXiv:1407.1974 (2014).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [accu,prediction] = kernel_knn_classification_new(test_kernel,train_label,K,test_label,options)
    ntest = size(test_kernel,2);
    if(size(test_kernel,1) ~= length(train_label))
        test_kernel = test_kernel';
        ntest = size(test_kernel,2);
    end
    test_label = test_label';
    train_label = train_label';
    nk = size(K,2);
    prediction = zeros(ntest,nk);
    for itest = 1:ntest
        cur_test = test_kernel(:,itest);
        [Y,I] = sort(cur_test,'descend');
        ttrain_label = train_label(I);
        for ik = 1:nk
            ktrain_label = ttrain_label(1:K(ik));

            ulabel = unique(ktrain_label);

            nlabel = size(ulabel,1);
            count = zeros(1,nlabel);
            for ilabel = 1:nlabel
                count(ilabel) = sum(ktrain_label==ulabel(ilabel));
            end
            [~,I] = max(count);
            prediction(itest,ik) = ulabel(I);
        end

       if options.verbose
            correct = (ulabel(I) == test_label(itest));
            fprintf('# Kernel kNN classifier: test:%03d, predict class: %03d --> ground truth :%03d (%d)\n', itest, ulabel(I), test_label(itest), correct);
       end

    end

    for ik = 1:nk
        accu(ik) = sum(prediction(:,ik) == test_label)/ntest;
    end
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function stokes_gnuplot ( header, n, x, y, u, v, s )

%*****************************************************************************80
%
%% STOKES_GNUPLOT writes the Stokes vector field to files for GNUPLOT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string HEADER, a header to be used to name the files.
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), Y(N), the coordinates of the evaluation points.
%
%    Input, real U(N), V(N), the velocity components.
%
%    Input, real S, a scale factor for the velocity vectors.
%

%
%  Write the data file.
%
  data_filename = strcat ( header, '_data.txt' );

  data_unit = fopen ( data_filename, 'wt' );

  for i = 1 : n
    fprintf ( data_unit, '  %10.4f  %10.4f  %10.4f  %10.4f\n', x(i), y(i), s * u(i), s * v(i) );
  end

  fclose ( data_unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data written to "%s".\n', data_filename );
%
%  Write the command file.
%
  command_filename = strcat ( header, '_commands.txt' );
  plot_filename = strcat ( header, '.png' );

  command_unit = fopen ( command_filename, 'wt' );

  fprintf ( command_unit, '#  %s\n', command_filename );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'set term png\n' );
  fprintf ( command_unit, 'set output "%s"\n', plot_filename );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, '#  Add titles and labels.\n' );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'set xlabel "<--- X --->"\n' );
  fprintf ( command_unit, 'set ylabel "<--- Y --->"\n' );
  fprintf ( command_unit, 'set title "Stokes velocity field"\n' );
  fprintf ( command_unit, 'unset key\n' );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, '#  Add grid lines.\n' );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'set grid\n' );
  fprintf ( command_unit, 'set size ratio -1\n' );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, '#  Timestamp the plot.\n' );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'set timestamp\n' );
  fprintf ( command_unit, 'plot "%s" using 1:2:3:4 with vectors \\\n', data_filename );
  fprintf ( command_unit, '  head filled lt 2 linecolor rgb "blue"\n' );
  fprintf ( command_unit, 'quit\n' );

  fclose ( command_unit );

  fprintf ( 1, '  Commands written to "%s".\n', command_filename );

  return
end


text:
function val = v_dtiRawPreprocessGE()
% Validate dti raw preprocessing of GE data
%
%    val = v_dtiRawPreprocessGE()
%
% This function checks that the average FA and B0 values resulting from
% dtiRawPreprocess are consistent with the expected value on different
% platforms (LINUX, WINDOWS, etc.).
%
% Requires vistadata and vistasoft on your path. The directory
% vistadata/diffusion/dtiRawPreprocess/GE contains all the necessary data
% files.
%
% This function does the following:
%
% (1) run dtiRawPreprocess on the files in the raw folder. This will align to the T1 and
% create a B0, BVECS and FA values (it creates more stuff but we only focus on these for the moment).
%
% (2) Show a montage of Alignment i.e., the T1 and DTI overalyed. Check for LR flips and correct alignment.
%
% (3) Compute FA, Mean diffusivity, radial diffusivity nd Axial diffusivity
%     across the brain and check the value obtained on different
%     platforms. THis will be done using: [fa,md,rd,ad] = dtiComputeFA(eigVal)
%
% (5) Load the stored FA, MD, RD, AD for the whole brain in:
%     GE/storedMeanDiffusionVals.mat
%
% (6) Recompute them from the data
%
% (7) Compute the difference between the stored and the recomputed ones.
%
% Example:
%  This function is meant to be called by mrvValidate, which will
% compute the difference between the stored values and the re-computed
% values.
%
%   mrvValidate([],[], 'v_dtiRawPreprocessGE');
%
% See also: mrvValidateAll.m
%
% FP and MP 7/6/2011
% Copyright Stanford team, mrVista, 2011

% This function takes too long to run. We want validation functions to be
% fast (say, less than 10 s). this one takes minutes to hours to run.
%
val = [];
return;

%% Get the data pathdata path
% Changing the last parameter and the function name changes
% scanner type (e.g., Siemens)
dataDir = fullfile(mrvDataRootPath,'diffusion','dtiRawPreprocess','GE');

%% Retain original directory, change to data directory
curDir = pwd;
cd(dataDir);

%% Run Preprocess
% setting clobber to 'always', so that output files will be silently replaced.
dtiRawPreprocess('raw/dti_g87_b1000.nii.gz', 't1.nii.gz',[],[],'always');

%% Show alignment to t1
% this was automatically computed by dtiRawPreprocess.m
imshow(imread('dti40trilin/t1pdd.png'));

%% Compute mean FA, RD, MD, AD values and check them with the stored one.

% load the stored values
% by convention mrvValidate.m loads up a storedVals data file
% with the same name of the validate function thta uses it
thisfunction = mfilename;
load(fullfile(mrvDataRootPath,'validate',[thisfunction(3:end),'.mat']));

% load the dti file.
dt = dtiLoadDt6('dti40trilin/dt6.mat');

% extract the eigen values.
eigVal = dt.dt6;

% compute the fractional, mean, radial and axial diffusivity
[fa,md,rd,ad] = dtiComputeFA(eigVal);

% compute the mean across the whole brain
val.fa = nanmean(fa(:));
val.md = nanmean(md(:));
val.rd = nanmean(rd(:));
val.ad = nanmean(ad(:));

% mrvValidate will use the output of this function to compare the
% recomputed values to the stored values.
% fr example doing something like this:
%
% check the computed with the stored values
% val.faErr = diff([mean_fa,meanVals.fa]);
% val.mdErr = diff([mean_md,meanVals.md]);
% val.rdErr = diff([mean_rd,meanVals.rd]);
% val.adErr = diff([mean_ad,meanVals.ad]);
%
% show results on matlab output
% errFields = fields(val);
% meanFields = fields(meanVals);
% for i = 1:length(fields(val))
%  fprintf('[%s] Error in ''%s'': %2.8f\n',mfilename, meanFields{i}, val.(errFields{i}));
% end

%% go back to the original directory, done!
cd(curDir)

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [Model, Info] = linear_sparse_stepwise(X,Y,Model,parm)
%  Estimate linear weight matrix for input-output mapping
%     Automatic Relevance Prior for each input dimension
%     is imposed to get sparse weight matrix
%
% Notice !!
%     Delay embedding is done in this module,
%     then input vector should be original input data without embedding
%
%   [Model, Info] = linear_sparse_stepwise(X,Y,Model,parm)
%     Scalar iteration version
% --- Input
%  Y  : Output data ( N x T x Ntrial )
%  X  : Input data  ( M x (T + (Dtau-1)*Tau)  x Ntrial)
%  N  =  # of output
%  M  =  # of input (original input space dimension)
%  T  =  # of time sample
%
%  Estimate the following model
%    Y(t) = W * [X(:, t + (Dtau-1)*Tau ); ...; X(:, t)]
%
%  Model : Structure for estimated model
%  if Model is empty, initialization is done before training
%  if Model is previous training result, re-training us done
%
%  parm  : Structure for learning parameter
%  parm.Npre_train :  # of VB-update in initial training
%  parm.Ntrain :  # of training
%  parm.Nskip  :  skip # for print
%  parm.a_min  :  Min value for pruning small variance component
%  parm.Prune  :  = 1 : Prune small variance & irrelevant input dimension
%
%  parm.Tau       = Lag time
%  parm.Dtau      = Number of embedding dimension
%    Total input dimension in embedding space is (M * parm.Dtau)
% --- Output
%  Model : Structure for estimated model
%  Model.SY  :  Noise variance         ( 1 x 1 )
%  Model.W   :  Weight matrix          ( N x M*D ) , D = parm.Dtau
%  Model.A   :  Prior weight variance  ( 1 x M*D )
%  Model.ix_act : Active index for W after pruning
%
%  Info  : Structure for learning process history
%  Info.FE  = LP + H : Free energy
%  Info.LP  = Log likelihood
%  Info.H   = - Model entropy
%
% 2009-10-18 Made by M. Sato

MINVAL = 1.0e-15;

fprintf('linear_sparse_stepwise start\n')

% Dimension
[N ,T  ,Ntrial]  = size(Y); % N = # of output
[M ,Tx ,Ntrialx] = size(X); % M = # of input without embedding

% Number of samples for each trial
if isfield(parm,'Trial'),
	Trial = parm.Trial;
	if length(Trial) ~= Ntrial, error('Trial sample seting is wrong');end
	if max(Trial) > T, error('Trial sample seting is wrong');end
else
	Trial = repmat(T, [Ntrial 1]);
end
% Set valid sample flag for each trial
trial_sw = zeros(N,T,Ntrial);
for n = 1:Ntrial
	trial_sw(:,1:Trial(n),n) = 1;
end

Tall = sum(Trial);

if Ntrial~=Ntrialx, error('# of trial is different in X and Y'); end;

Nskip  = 100;   % skip steps for display info
a_min = 1e-14;	% Minimum value for weight pruning
Fdiff = 1e-12; % Threshold for convergence
Ncheck = 100;   % Minimum number of training iteration
Fstep  = 5;     % Free energy convergence check step
Prune  = 1;     % Prune mode
space_ARD=1;	% if space_ARD=1, ARD is done only for space dimension

if isfield(parm,'Nskip'), Nskip  = parm.Nskip; end;
if isfield(parm,'Fdiff'), Fdiff   = parm.Fdiff; end;
if isfield(parm,'a_min'), a_min   = parm.a_min ; end;
if isfield(parm,'Prune'), Prune = parm.Prune; end;
if isfield(parm,'Ncheck'), Ncheck = parm.Ncheck; end
if isfield(parm,'Fstep'),Fstep  = parm.Fstep ; end
if isfield(parm,'Fstep'),Fstep  = parm.Fstep ; end
if isfield(parm,'space_ARD'),space_ARD  = parm.space_ARD ; end

% # of embedding dimension
if isfield(parm,'Dtau')
	D    = parm.Dtau;
	tau  = parm.Tau;
else
	D    = 1;
	tau  = 1;
end
% Parameters
Ntrain = parm.Ntrain;

if Ntrain < 1, Info = []; return; end;

if isfield(parm,'Npre_train')
	Npre_train = parm.Npre_train;
else
%	Npre_train = Ntrain;
	if Tall >= 2*M*D
		Npre_train = 0;
	elseif Tall >= M*D
		Npre_train = fix(Ntrain/2);
	else
		Npre_train = Ntrain;
	end
end

if Npre_train > Ntrain, Npre_train = Ntrain; end;

if isfield(parm,'Nupdate')
	Nupdate = parm.Nupdate;
else
	Nupdate = 1;
end

fprintf('--- Output Dimension  = %d\n',N)
fprintf('--- Input  Dimension  = %d\n',M)
fprintf('--- Embedding  Dimension  = %d\n',D)
fprintf('--- Number of trials  = %d\n',Ntrial)
fprintf('--- Number of training sample = %d\n',Tall)
fprintf('--- Total update iteration    = %d (%d)\n',Ntrain,Npre_train)

% original input dimension
Xdim  = M;
M = Xdim*D;
M_ALL = M;

%
% --- Initialization
%

% Input/Output variance
sx = mean((X(:) - mean(X(:))).^2);
sy = mean((Y(:) - mean(Y(:))).^2);

SY0 = mean(sy);
A0  = 1./mean(sx);

if isfield(parm,'Ta0') && parm.Ta0 > 0,
	Ta0 = parm.Ta0;
	a0  = parm.a0 * A0;
else
	Ta0 = 0;
	a0  = 1;
end

% Input covariance
Trial = Trial+(D-1)*tau;

XX = zeros(Xdim,1);
for n=1:Ntrial
	XX = XX + sum(X(:,1:Trial(n), n).^2 , 2);
end
XX  = XX / sum(Trial);
%XX  = sum(sum(X.^2,3),2)/(Tx*Ntrial);
XX  = repmat(XX', [1 D]);% 1 x M

if isfield(Model,'ix_act')
	SY  = mean(Model.SY);  % 1 x 1

	% Active index
	ix_act = Model.ix_act;
	A = zeros(1,M_ALL);
	W = zeros(N,M_ALL);

	A(ix_act)   = sum(Model.A,1);
	W(:,ix_act) = Model.W;

	if M_ALL ~= Xdim*D,
		fprintf('M_ALL=%d,Xdim=%d,D=%d\n',M_ALL ,Xdim ,D)
		error('M_ALL ~= Xdim*D')
	end

	% Active index in input space without embedding
	IX_dim = find( sum(reshape(A,[Xdim,D]),2) > 0);

	IX_act = repmat( (0:(D-1))* Xdim ,[length(IX_dim) 1]) ...
	       + repmat(IX_dim, [1 D]);
	IX_act = IX_act(:);

	W  = W(:,IX_act) ;  % N x (M*D)
	A  = A(IX_act) ;  % 1 x (M*D)
	X  = X(IX_dim,:,:);
	XX = XX(IX_act);

	M     = length(IX_act);
	Xdim  = length(IX_dim);
else
	A = repmat(A0, [1,M_ALL]);
	W = zeros(N,M_ALL);
	SY  = SY0;
	% Save original input for pruning
	IX_act  = 1:M_ALL;
	IX_dim  = 1:Xdim;
end

ix_act  = 1:M;
ix_dim  = 1:Xdim;
M_all   = M;

% Delay embedding index for W
Wid = [(0:(D-1))'* Xdim + 1 , (1:D)'* Xdim];

% Working variable
G_A = zeros(N,M);       % N x M

if size(A,1)==Xdim && size(A,2)==1
	A   = repmat(A',[1 D]);
elseif size(A,1)==1 && size(A,2)==Xdim
	A   = repmat(A,[1 D]);
end

SW = SY./( Tall*XX + SY./A );

fprintf('a_min = %g\n', a_min)
fprintf('SY0 = %g\n', SY0)
fprintf('SY  = %g\n', SY)

% Free energy histry
FE  = zeros(Ntrain,1);
LP  = zeros(Ntrain,1);
H   = zeros(Ntrain,1);
Err = zeros(Ntrain,1);
Mhist = zeros(Ntrain,1);

% ARD hyper param. history
if isfield(parm,'Debug') && ~isempty(parm.Debug) && parm.Debug > 0
	Debug = 1;
	A_tmp = zeros(N*M_ALL, ceil(Ntrain/Nskip));
else
	Debug = 0;
end

% recover all component
A_all  = zeros(1,M_all);       % N x M

dY = error_delay_time_sw(X,Y,W,tau,Trial);
dY = dY .* trial_sw;                           % mask invalid samples

k_save  = 0;

%%%%%% Learning Loop %%%%%%
for k=1:Ntrain
	% Ainv = alpha , A = 1/alpha
	% (Y-W*X)^2/SY + W^2/A
	% ( (Y-W*X)^2 + W^2*Ainv )/SY
	Ainv    = SY./A;

    W   = weight_update_embed_sw(X, dY, W, Tall*XX, Ainv, tau, Trial);
	dY  = error_delay_time_sw(X, Y, W, tau, Trial);
	dY  = dY .* trial_sw;                           % mask invalid samples

    dYY = sum(dY(:).^2)/(N*Tall);
	WW  = sum(W.^2,1);

    % Noise variance update
    SY  = dYY + sum(WW .* Ainv)/(N*Tall);
    % Prevent zero variance
    SY  = max( SY, MINVAL);

    % Log variance
    SWA     = max( SW .* Ainv , MINVAL);
    log_sw  = N*(sum( log(SWA) - SWA + 1 ));
    log_sy  = N*( log(SY) );
    log_a   = Ta0*(sum(log(Ainv) - a0.*Ainv + 1));

    % Free energy
    LP(k)  = - (0.5*Tall) * log_sy ;
    H(k)   = 0.5*( log_sw + log_a - sum(WW .*Ainv) );
    FE(k)  = LP(k) + H(k);
    Err(k) = sum(dYY)/(N*SY0);

    % Weight variance

    SW = SY./( Tall*XX + Ainv );

    % Hyper parameter for weight variance (ARD)
	if mod(k,Nupdate) == 0,
	    % A = Alpha^-1 * SY
		% Ainv = SY./A;
	    % SW = SY./( Tall*XX + Ainv );

	    % G_A = 1 - SW./A;
	    %     = (Tall*XX + Ainv - SY/A) * SW/SY
	    G_A = Tall.*XX.*SW./SY;
	    G_A = max((G_A), MINVAL);

		% N*A = N * (Alpha * SY) =  (W.^2) + N * SW  ;
		if k <= Npre_train,
			% VB update rule (Stable)
			%  ARD for each weight
			A  = (WW + N*SW + 2*Ta0*a0)/( N + 2*Ta0 );
		else
			% Accelerated update rule (Unstable)
		    A  = sqrt( A.* WW./(G_A * N));
		    %A  = (WW + 2*Ta0*a0)./(G_A * N + 2*Ta0);
		end

		if space_ARD==1
			Am = mean(reshape(A,[M/D,D]),2);
			A  = repmat(Am', 1,D);
		end

	    % Prune small variance
	    if Prune == 1
		    % Find active input dimension
		    ix_act_old = ix_act;
		    ix_dim_old = ix_dim;

		    % recover all component
			A_all  = zeros(1,M_all);       % 1 x M
		    %  A_all(:,ix_act)  = A;
		    A_all(ix_act) = WW/max(WW);

		    % find active input dimension (absolute index)
		    ix_dim = find( sum(reshape(A_all,[Xdim,D]),2) > a_min );
		    ix_act = repmat(ix_dim, [1 D]) ...
		           + repmat([0:D-1]*Xdim, [length(ix_dim) 1]);
		    ix_act = ix_act(:);

		    Mnew   = length(ix_act);  		% # of effective input

		    if Mnew < M,
			    % convert to relative index
			    jx_act = trans_index(ix_act,ix_act_old,M_all);
			    jx_dim = trans_index(ix_dim,ix_dim_old,Xdim);

			    M   = Mnew;
			    A   = A(jx_act) ;  % N x M
			    W   = W(:,jx_act) ;  % N x M
			    SW  = SW(jx_act);  % N x M

			    X	 = X(jx_dim,:,:);  	 	% M x T
			    XX	 = XX(jx_act);  	 	% M x T
			end
	    end
	    % END of if Prune == 1

		A = max(A,MINVAL);

	end
	% END of if mod(k,Nupdate) == 0

	Mhist(k) = M;
	SW  = SW /SY;

    if mod(k, Nskip)==0
        % Save history
		if Debug == 1
        	k_save = k_save + 1;
        	A_tmp(:,k_save) = A_all(:);
		end

        fprintf('Iter = %4d, M = %4d, err = %g, F = %g, SY= %g, H = %g\n', ...
               k, length(ix_act), Err(k), FE(k), SY, - 2*H(k)/log(T));
    end

    % Convergence check
	if k > Ncheck,
		Fdif = (FE(k) - FE(k-Fstep))/(abs(FE(k))+eps);
	else
		Fdif = Fdiff + 1;
	end

	if (Fdiff > abs(Fdif)),
		fprintf('Converged : Free energy change = %g\n',Fdif)
		break;
	end;
end

ix_act = IX_act(ix_act);

% Active index
Model.ix_act = ix_act;
Model.M_all  = M_ALL ;

% Save output variable
Model.A    = A ;
Model.W    = W ;
Model.SY   = SY;
Model.SW   = SW; % = 1./(Tall*XX + diag(Ainv))

Model.method = 'linear_sparse_stepwise';
Model.mode   = 'scalar';
Model.sparse = 'sparse';

% Save history
Info.FE  = FE(1:k);
Info.LP  = LP(1:k);
Info.H   = H(1:k) ;
Info.Err = Err(1:k);
Info.M   = Mhist(1:k);

if exist('A_tmp','var')
	Info.A   = A_tmp(:,1:k_save) ;
end

% recover all component
%W_all   = zeros(N,M_all);
%SW_all  = zeros(N,M_all);
%
%A_all(:,ix_act)  = A;
%W_all(:,ix_act)  = W ;  % N x M
%SW_all(:,ix_act) = SW;  % N x M

%%% ---- Index transformation from old active_index to current active_index
function	jx = trans_index(ix,ix_old,M)

N = length(ix_old);
Itrans = zeros(M,1);
Itrans(ix_old) = 1:N;

jx = Itrans(ix);

author: KamitaniLab

repo: GenericObjectDecoding


text:
function [x ft] = EProjSimplex_new(v, k)

%
%% Problem
%
%  min  1/2 || x - v||^2
%  s.t. x>=0, 1'x=1
%

if nargin < 2
    k = 1;
end;

ft=1;
n = length(v);

v0 = v-mean(v) + k/n;
%vmax = max(v0);
vmin = min(v0);
if vmin < 0
    f = 1;
    lambda_m = 0;
    while abs(f) > 10^-10
        v1 = v0 - lambda_m;
        posidx = v1>0;
        npos = sum(posidx);
        g = -npos;
        f = sum(v1(posidx)) - k;
        lambda_m = lambda_m - f/g;
        ft=ft+1;
        if ft > 100
            x = max(v1,0);
            break;
        end;
    end;
    x = max(v1,0);
else
    x = v0;
end;

author: CHLWR


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% TLD is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with TLD.  If not, see <http://www.gnu.org/licenses/>.

function tld = tldLearning(tld,I)

bb    = tld.bb(:,I); % current bounding box
img   = tld.img{I}; % current image

% Check consistency -------------------------------------------------------

pPatt  = tldGetPattern(img,bb,tld.model.patchsize); % get current patch
[pConf1,~,pIsin] = tldNN(pPatt,tld); % measure similarity to model

if pConf1 < 0.5, disp('Fast change.'); tld.valid(I) = 0; return; end % too fast change of appearance
if var(pPatt) < tld.var, disp('Low variance.'); tld.valid(I) = 0; return; end % too low variance of the patch
if pIsin(3) == 1, disp('In negative data.'); tld.valid(I) = 0; return; end % patch is in negative data

% Update ------------------------------------------------------------------

% generate positive data
overlap  = bb_overlap(bb,tld.grid); % measure overlap of the current bounding box with the bounding boxes on the grid
[pX,pEx] = tldGeneratePositiveData(tld,overlap,img,tld.p_par_update); % generate positive examples from all bounding boxes that are highly overlappipng with current bounding box
pY       = ones(1,size(pX,2)); % labels of the positive patches

% generate negative data
idx      = overlap < tld.n_par.overlap & tld.tmp.conf >= 1; % get indexes of negative bounding boxes on the grid (bounding boxes on the grid that are far from current bounding box and which confidence was larger than 0)
overlap  = bb_overlap(bb,tld.dt{I}.bb); % measure overlap of the current bounding box with detections
nEx      = tld.dt{I}.patch(:,overlap < tld.n_par.overlap); % get negative patches that are far from current bounding box

fern(2,[pX tld.tmp.patt(:,idx)],[pY zeros(1,sum(idx))],tld.model.thr_fern,2); % update the Ensemble Classifier (reuses the computation made by detector)
tld = tldTrainNN(pEx,nEx,tld); % update nearest neighbour

repo: OpenTLD


text:
function [logp, yhat, res] = tapas_beta_obs(r, infStates, ptrans)
% Calculates the log-probability of responses representing probabilities on the unit interval
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013-2016 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Transform nu-prime to its native space
nupr = exp(ptrans(1));

% Initialize returned log-probabilities, predictions,
% and residuals as NaNs so that NaN is returned for all
% irregualar trials
n = size(infStates,1);
logp = NaN(n,1);
yhat = NaN(n,1);
res  = NaN(n,1);

% Predictions or posteriors?
if strcmp(r.c_prc.model,'tapas_rw_binary')
    mu = tapas_sgm(infStates(:,1), 1);
else
    mu = infStates(:,1,1); % Default: predictions (ie, mu1hat)
    if r.c_obs.predorpost == 2
        mu = tapas_sgm(infStates(:,2,3), 1); % Alternative: posteriors (ie, sgm(mu2))
    end
end

% Special cases
if strcmp(r.c_prc.model,'hgf_whichworld')
    mu = tapas_sgm(infStates(:,2,1,3), 1);
end
if strcmp(r.c_prc.model,'ph_binary')
    mu = infStates(:,2);
end

% Weed irregular trials out from inferred states and responses
mu(r.irr) = [];
y = r.y(:,1);
y(r.irr) = [];

% y has to be in the *open* unit interval
%y(y==0) = 1e-4;
%y(y==1) = 1-1e-4;
y = 0.95.*(y-0.5)+0.5; % Shrink all y values toward 1/2 by a factor of 0.95

% Nu is nu-prime plus two (sometimes)
%nu = nupr+2;
nu = nupr;

% Calculate alpha and beta from mu and nu
al = mu.*nu;
be = nu - al;

% Calculate log-probabilities for non-irregular trials
reg = ~ismember(1:n,r.irr);
logp(reg) = log(betaDens(y,al,be));
yhat(reg) = mu;
res(reg) = y-mu;

end

function p = betaDens(x,alpha,beta)
% Check whether x is in the unit interval
if any(x(:)<0) || any(x(:)>1)
    error('tapas:hgf:BetaObs:ArgNotInUnitIntrv', 'Error: first argument to betaDens must be in the unit interval.');
end
% Check whether alpha and beta are greater than 0
if any(alpha(:)<0) || any(beta(:)<0)
    error('tapas:hgf:BetaObs:AlphaOrBetaNeg', 'Error: alpha and beta have to be non-negative.');
end
% Calculate beta density
p = gamma(alpha+beta)./(gamma(alpha).*gamma(beta)).*x.^(alpha-1).*(1-x).^(beta-1);
end

author: translationalneuromodeling

repo: tapas