diff --git a/lib/Matrix.pm b/lib/Matrix.pm index 18c8093ff9..af6bc8073d 100644 --- a/lib/Matrix.pm +++ b/lib/Matrix.pm @@ -1,11 +1,17 @@ =head1 NAME -Matrix - Matrix of Reals +lib/Matrix - Matrix of Reals -Implements overrides for MatrixReal.pm for WeBWorK =head1 DESCRIPTION +Implements overrides for MatrixReal.pm for WeBWorK +In general it is better to use MathObjects Matrices (Value::Matrix) +in writing PG problem. The answer checking is much superior with better +error messages for syntax errors in student entries. Some of the +subroutines in this file are still used behind the scenes +by Value::Matrix to perform calculations, +such as decompose_LR(). =head1 SYNOPSIS @@ -68,8 +74,22 @@ sub _stringify { return($s); } +=head3 Accessor functions + + (these are deprecated for direct use. Use the covering Methods + provided by MathObject Matrices instead.) + + L($matrix) - return matrix L of the LR decomposition + R($matrix) - return matrix R of the LR decomposition + PL($matrix) return + PR($matrix + +Original matrix is P_L * L * R *P_R # obtain the Left Right matrices of the decomposition and the two pivot permutation matrices # the original is M = PL*L*R*PR + +=cut + sub L { my $matrix = shift; my $rows = $matrix->[1]; @@ -119,10 +139,14 @@ sub PR { # use this permuation on the right PL*L*R*PR =M } # obtain the Left Right matrices of the decomposition and the two pivot permutation matrices # the original is M = PL*L*R*PR -=head4 + + +=item rh_options Method $matrix->rh_options +Meant for internal use when dealing with MatrixReal1 + =cut sub rh_options { @@ -132,14 +156,18 @@ sub rh_options { $self->[$MatrixReal1::OPTION_ENTRY]; # provides a reference to the options hash MEG } -=head4 - Method $matrix->trace - +=item trace + + Method: $matrix->trace Returns: scalar which is the trace of the matrix. + +Used by MathObject Matrices for calculating the trace. +Deprecated for direct use in PG questions. =cut + sub trace { my $self = shift; my $rows = $self->[1]; @@ -152,9 +180,12 @@ sub trace { $sum; } -=head4 - Method $matrix->new_from_array_ref +=item new_from_array_ref + + Method $new_matrix = $matrix->new_from_array_ref ([[a,b,c],[d,e,f]]) + +Deprecated in favor of using creation tools for MathObject Matrices =cut @@ -168,10 +199,12 @@ sub new_from_array_ref { # this will build a matrix or a row vector from [a, b $matrix; } -=head4 +=item array_ref Method $matrix->array_ref +Converts Matrix from an ARRAY to an ARRAY reference. + =cut sub array_ref { @@ -179,10 +212,12 @@ sub array_ref { $this->[0]; } -=head4 +=item list Method $matrix->list +Converts a Matrix column vector to an ARRAY (list). + =cut sub list { # this is used only for column vectors @@ -196,29 +231,14 @@ sub list { # this is used only for column vectors @list; } -=head4 - - Method $matrix->new_from_list - -=cut - -sub new_from_list { # this builds a row vector from an array - my $class = shift; - my @list = @_; - my $cols = @list; - my $rows = 1; - my $matrix = new Matrix($rows, $cols); - my $i=1; - while(@list) { - my $elem = shift(@list); - $matrix->assign($i++,1, $elem); - } - $matrix; -} -=head4 +=item new_row_matrix Method $matrix->new_row_matrix + +Deprecated -- there are better tools for MathObject Matrices. + +Create a row 1 by n matrix from a list. This subroutine appears to be broken =cut @@ -236,10 +256,12 @@ sub new_row_matrix { # this builds a row vector from an array $matrix; } -=head4 +=item proj Method $matrix->proj - + Provides behind the scenes calculations for MathObject Matrix->proj + Deprecated for direct use in favor of methods of MathObject matrix + =cut sub proj{ @@ -248,9 +270,12 @@ sub proj{ $self * $self ->proj_coeff($vec); } -=head4 +=item proj_coeff Method $matrix->proj_coeff + +Provides behind the scenes calculations for MathObject Matrix->proj_coeff +Deprecated for direct use in favor of methods of MathObject matrix =cut @@ -267,10 +292,12 @@ sub proj_coeff{ $x_vector; } -=head4 +=item new_column_matrix Method $matrix->new_column_matrix +Create column matrix from an ARRAY reference (list reference) + =cut sub new_column_matrix { @@ -286,13 +313,15 @@ sub new_column_matrix { $matrix; } -=head4 - - This method takes an array of column vectors, or an array of arrays, - and converts them to a matrix where each column is one of the previous - vectors. +=item new_from_col_vecs Method $matrix->new_from_col_vecs + +Deprecated: The tools for creating MathObjects Matrices are simpler. +This method takes an array of column vectors, or an array of arrays, +and converts them to a matrix where each column is one of the previous +vectors. + =cut @@ -341,9 +370,11 @@ sub new_from_col_vecs =cut -=head4 +=item cp - Method $matrix->new_from_col_vecs + Function: cp() + +Provides ability to use perl complex numbers. N =cut @@ -354,7 +385,7 @@ sub cp { # MEG makes new copies of complex number return $w; } -=head4 +=item copy Method $matrix->copy @@ -397,7 +428,7 @@ sub copy # MEG added 6/25/03 to accomodate complex entries -=head4 +=item conj Method $matrix->conj @@ -409,7 +440,7 @@ sub conj { $elem; } -=head4 +=item transpose Method $matrix->transpose @@ -458,10 +489,13 @@ sub transpose $matrix1; } -=head4 +=item decompose_LR Method $matrix->decompose_LR +Used by MathObjects Matrix for LR decomposition +Deprecated for direct use in PG problems. + =cut sub decompose_LR diff --git a/lib/Value/Matrix.pm b/lib/Value/Matrix.pm index e853ff5ff5..1cbbdcd249 100644 --- a/lib/Value/Matrix.pm +++ b/lib/Value/Matrix.pm @@ -3,6 +3,111 @@ # Implements the Matrix class. # # @@@ Still needs lots of work @@@ + +=head1 NAME + + Value::Matrix class + + +References: + +MathObject Matrix methods: L +MathObject Contexts: L +CPAN RealMatrix docs: L + +Allowing Matrices in Fractions: +L + + Context()->parens->set("[" => {formMatrix => 1}); + +Files interacting with Matrices: + +L L +L + +L -- checking whether vectors form a basis +L -- tools for row reduction via elementary matrices +L -- Generates unimodular matrices with real entries +L +L +L +L +L +L +Contexts + + Matrix -- allows students to enter [[3,4],[3,6]] + -- formMatrix =>1 also allows this? + Complex-Matrix -- allows complex entries + +Creation methods + + $M1 = Matrix([1,2],[3,4]); + $M2 = Matrix([5,6],[7,8]); + $v = Vector(9,10); + $w = ColumnVector(9,10); # differs in how it is printed + +Commands added in Value::matrix + + Conversion + $matrix->values produces [[3,4,5],[1,3,4]] recursive array references of numbers (not MathObjects) + $matrix->wwMatrix produces CPAN MatrixReal1 matrix, used for computation subroutines + + Information + $matrix->dimension: ARRAY + + Access values + + row : MathObjectMatrix + column : MathObjectMatrix + element : Real or Complex value + + Assign values + + these need to be added: + +see C in MatrixReduce and L + + Advanced + $matrix->data: ARRAY reference (internal data) of MathObjects (Real,Complex, Fractions) + stored at each location. + + +Passthrough methods covering subroutines in Matrix.pm which overrides or +augment CPAN's MatrixReal1.pm. Matrix is a specialized subclass of MatrixReal1.pm + +The actual calculations are done in Matrix.pm + + trace + proj + proj_coeff + L + R + PL + PR + +Passthrough methods covering subroutines in MatrixReal1.pm +(this has been modified to handle complex numbers) +The actual calculations are done in MatrixReal1.pm subroutines + + condition + det + inverse + is_symmetric + decompose_LR + dim + norm_one + norm_max + kleene + normalize + solve_LR (also solve()) + order_LR (also order() + solve_GSM + solve_SSM + solve_RM + +=cut + # package Value::Matrix; my $pkg = 'Value::Matrix'; @@ -17,7 +122,7 @@ our @ISA = qw(Value); # a point, vector or matrix object, a matrix-valued formula, or a string # that evaluates to a matrix # -sub new { +sub new { #internal my $self = shift; my $class = ref($self) || $self; my $context = (Value::isContext($_[0]) ? shift : $self->context); my $M = shift; $M = [] unless defined $M; $M = [$M,@_] if scalar(@_) > 0; @@ -38,7 +143,7 @@ sub new { # (Recursively) make a matrix from a list of array refs # and report errors about the entry types # -sub matrixMatrix { +sub matrixMatrix { #internal my $self = shift; my $class = ref($self) || $self; my $context = shift; my ($x,$m); my @M = (); my $isFormula = 0; @@ -62,7 +167,7 @@ sub matrixMatrix { # Form a 1 x n matrix from a list of numbers # (could become a row of an m x n matrix) # -sub numberMatrix { +sub numberMatrix { #internal my $self = shift; my $class = ref($self) || $self; my $context = shift; my @M = (); my $isFormula = 0; @@ -209,7 +314,7 @@ sub mult { # # Constant multiplication # - if (Value::matchNumber($r) || Value::isComplex($r)) { + if (Value::isNumber($r)) { my @coords = (); foreach my $x (@{$l->data}) {push(@coords,$x*$r)} return $self->make(@coords); @@ -247,8 +352,7 @@ sub mult { sub div { my ($l,$r,$flag) = @_; my $self = $l; Value::Error("Can't divide by a Matrix") if $flag; - Value::Error("Matrices can only be divided by Numbers") - unless (Value::matchNumber($r) || Value::isComplex($r)); + Value::Error("Matrices can only be divided by Numbers") unless Value::isNumber($r); Value::Error("Division by zero") if $r == 0; my @coords = (); foreach my $x (@{$l->data}) {push(@coords,$x/$r)} diff --git a/lib/Value/Point.pm b/lib/Value/Point.pm index 6ab10d25bf..9d508ef265 100644 --- a/lib/Value/Point.pm +++ b/lib/Value/Point.pm @@ -81,8 +81,7 @@ sub sub { sub mult { my ($l,$r) = @_; my $self = $l; - Value::Error("Points can only be multiplied by Numbers") - unless (Value::matchNumber($r) || Value::isComplex($r)); + Value::Error("Points can only be multiplied by Numbers") unless Value::isNumber($r); my @coords = (); foreach my $x ($l->value) {push(@coords,$x*$r)} return $self->make(@coords); @@ -91,8 +90,7 @@ sub mult { sub div { my ($l,$r,$flag) = @_; my $self = $l; Value::Error("Can't divide by a Point") if $flag; - Value::Error("Points can only be divided by Numbers") - unless (Value::matchNumber($r) || Value::isComplex($r)); + Value::Error("Points can only be divided by Numbers") unless Value::isNumber($r); Value::Error("Division by zero") if $r == 0; my @coords = (); foreach my $x ($l->value) {push(@coords,$x/$r)} diff --git a/lib/Value/Vector.pm b/lib/Value/Vector.pm index 16116a2dab..52ac31f79f 100644 --- a/lib/Value/Vector.pm +++ b/lib/Value/Vector.pm @@ -92,8 +92,7 @@ sub sub { sub mult { my ($l,$r,$flag) = @_; my $self = $l; - Value::Error("Vectors can only be multiplied by Numbers") - unless (Value::matchNumber($r) || Value::isComplex($r)); + Value::Error("Vectors can only be multiplied by Numbers") unless Value::isNumber($r); my @coords = (); foreach my $x ($l->value) {push(@coords,$x*$r)} return $self->make(@coords); @@ -102,8 +101,7 @@ sub mult { sub div { my ($l,$r,$flag) = @_; my $self = $l; Value::Error("Can't divide by a Vector") if $flag; - Value::Error("Vectors can only be divided by Numbers") - unless (Value::matchNumber($r) || Value::isComplex($r)); + Value::Error("Vectors can only be divided by Numbers") unless Value::isNumber($r); Value::Error("Division by zero") if $r == 0; my @coords = (); foreach my $x ($l->value) {push(@coords,$x/$r)} diff --git a/macros/MatrixCheckers.pl b/macros/MatrixCheckers.pl index 6f4c963173..ec13808365 100644 --- a/macros/MatrixCheckers.pl +++ b/macros/MatrixCheckers.pl @@ -87,23 +87,23 @@ =head1 DESCRIPTION are produced by C<\(\Bigg\lbrace\)> and C<\(\Bigg\rbrace\)>, are a matter of personal preference (since a basis is an ordered set, I like to include braces). -=over 12 -Context()->texStrings; -BEGIN_TEXT -Find an orthonormal basis for... -$BR -$BR -$BCENTER -\(\Bigg\lbrace\) -\{ $multians->ans_array(15) \}, -\{ $multians->ans_array(15) \} -\(\Bigg\rbrace.\) -$ECENTER -END_TEXT -Context()->normalStrings; -=back + Context()->texStrings; + BEGIN_TEXT + Find an orthonormal basis for... + $BR + $BR + $BCENTER + \(\Bigg\lbrace\) + \{ $multians->ans_array(15) \}, + \{ $multians->ans_array(15) \} + \(\Bigg\rbrace.\) + $ECENTER + END_TEXT + Context()->normalStrings; + + The answer evaluation section of the PG file is totally standard. diff --git a/macros/MatrixReduce.pl b/macros/MatrixReduce.pl index 010d732269..a77925dc76 100644 --- a/macros/MatrixReduce.pl +++ b/macros/MatrixReduce.pl @@ -14,27 +14,45 @@ =head1 SYNOPSIS =over 12 -=item Get the reduced row echelon form: C<$Areduced = rref($A);> Should be used in the fraction context with all entries of $A made into fractions. +=item Get the reduced row echelon form: C<$Areduced = rref($A);> + +Should be used in the fraction context with all entries of $A made into fractions. -=item Make matrix entries do fraction arithmetic (rather than decimal arithmetic): After selecting the Fraction context using Context('Fraction')->parens->set("[" => {formMatrix => 1}), C<$A = apply_fraction_to_matrix_entries($A);> applies Fraction() to all of the entries of $A, which makes subsequent matrix algebra computations with $A use fraction arithmetic. +=item Make matrix entries do fraction arithmetic (rather than decimal arithmetic): + +After selecting the Fraction context using Context('Fraction')->parens->set("[" => {formMatrix => 1}), C<$A = apply_fraction_to_matrix_entries($A);> applies Fraction() to all of the entries of $A, which makes subsequent matrix algebra computations with $A use fraction arithmetic. =item Get the reduced column echelon form: C<$Areduced = rcef($A);> -=item Change the value of a matrix entry: C changes the [2,3] entry to the value 50. +=item Change the value of a matrix entry: C + +changes the [2,3] entry to the value 50. =item Construct an n x n identity matrix: C<$E = identity_matrix(5);> -=item Construct an n x n elementary matrix that will permute rows i and j: C<$E = elem_matrix_row_switch(5,2,4);> creates a 5 x 5 identity matrix and swaps rows 2 and 4. +=item Construct an n x n elementary matrix that will permute rows i and j: + +C<$E = elem_matrix_row_switch(5,2,4);> creates a 5 x 5 identity matrix and swaps rows 2 and 4. + +=item Construct an n x n elementary matrix that will multiply row i by s: C<$E = elem_matrix_row_mult(5,2,4);> + +creates a 5 x 5 identity matrix and swaps puts 4 in the second spot on the diagonal. -=item Construct an n x n elementary matrix that will multiply row i by s: C<$E = elem_matrix_row_mult(5,2,4);> creates a 5 x 5 identity matrix and swaps puts 4 in the second spot on the diagonal. +=item Construct an n x n elementary matrix that will multiply row i by s: C<$E3 = elem_matrix_row_add(5,3,1,35);> -=item Construct an n x n elementary matrix that will multiply row i by s: C<$E3 = elem_matrix_row_add(5,3,1,35);> creates a 5 x 5 identity matrix and swaps puts 35 in the (3,1) position. +creates a 5 x 5 identity matrix and swaps puts 35 in the (3,1) position. -=item Perform the row switch transform that swaps (row i) with (row j): C<$Areduced = row_switch($A,2,4);> swaps rows 2 and 4 in matrix $A. +=item Perform the row switch transform that swaps (row i) with (row j): C<$Areduced = row_switch($A,2,4);> -=item Perform the row multiplication transform s * (row i) placed into (row i): C<$Areduced = row_mult(A,2,10);> multiplies every entry in row 2 of $A by 10. +swaps rows 2 and 4 in matrix $A. -=item Perform the row addition transform (row i) + s * (row j) placed into (row i): C<$Areduced = row_add($A,2,1,10);> adds 10 times row 1 to row 2 and places the result in row 2. (Same as constructing $E to be the identity with 10 placed in entry (2,1), then multiplying $E * $A.) +=item Perform the row multiplication transform s * (row i) placed into (row i): C<$Areduced = row_mult(A,2,10);> + +multiplies every entry in row 2 of $A by 10. + +=item Perform the row addition transform (row i) + s * (row j) placed into (row i): C<$Areduced = row_add($A,2,1,10);> + +adds 10 times row 1 to row 2 and places the result in row 2. (Same as constructing $E to be the identity with 10 placed in entry (2,1), then multiplying $E * $A.) =back @@ -42,61 +60,59 @@ =head1 DESCRIPTION Usage: -=over 12 - -DOCUMENT(); -loadMacros( -"PGstandard.pl", -"MathObjects.pl", -"MatrixReduce.pl", # automatically loads contextFraction.pl and MathObjects.pl -"PGcourse.pl", -); -$showPartialCorrectAnswers = 0; -TEXT(beginproblem()); - -# Context('Matrix'); # for decimal arithmetic -Context('Fraction'); # for fraction arithmetic - -$A = Matrix([ -[random(-5,5,1),random(-5,5,1),random(-5,5,1),3], -[random(-5,5,1),random(-5,5,1),random(-5,5,1),0.75], -[random(-5,5,1),random(-5,5,1),random(-5,5,1),9/4], -]); - -$A = apply_fraction_to_matrix_entries($A); # try commenting this line out for different results - -$Arref = rref($A); - -$Aswitch = row_switch($A, 2, 3); - -$Amult = row_mult($A, 2, 4); - -$Aadd = row_add($A, 2, 1, 10); - -$E = elem_matrix_row_add(3,2,1,10); -$EA = $E * $A; - -$E1 = elem_matrix_row_switch(5,2,4); -$E2 = elem_matrix_row_mult(5,4,Fraction(1/10)); -$E3 = elem_matrix_row_add(5,3,1,35); -$E4 = identity_matrix(4); -change_matrix_entry($E4,[3,2],10); - -Context()->texStrings; -BEGIN_TEXT -The original matrix and its row reduced echelon form: -\[ $A \sim $Arref. \] -$BR -The original matrix with rows switched, multiplied, or added together: -\[ $Aswitch, $Amult, $Aadd. \] -$BR -Some elementary matrices. -\[$E1, $E2, $E3, $E4\] -END_TEXT -Context()->normalStrings; - -COMMENT('MathObject version.'); -ENDDOCUMENT(); + DOCUMENT(); + loadMacros( + "PGstandard.pl", + "MathObjects.pl", + "MatrixReduce.pl", # automatically loads contextFraction.pl and MathObjects.pl + "PGcourse.pl", + ); + $showPartialCorrectAnswers = 0; + TEXT(beginproblem()); + + # Context('Matrix'); # for decimal arithmetic + Context('Fraction'); # for fraction arithmetic + + $A = Matrix([ + [random(-5,5,1),random(-5,5,1),random(-5,5,1),3], + [random(-5,5,1),random(-5,5,1),random(-5,5,1),0.75], + [random(-5,5,1),random(-5,5,1),random(-5,5,1),9/4], + ]); + + $A = apply_fraction_to_matrix_entries($A); # try commenting this line out for different results + + $Arref = rref($A); + + $Aswitch = row_switch($A, 2, 3); + + $Amult = row_mult($A, 2, 4); + + $Aadd = row_add($A, 2, 1, 10); + + $E = elem_matrix_row_add(3,2,1,10); + $EA = $E * $A; + + $E1 = elem_matrix_row_switch(5,2,4); + $E2 = elem_matrix_row_mult(5,4,Fraction(1/10)); + $E3 = elem_matrix_row_add(5,3,1,35); + $E4 = identity_matrix(4); + change_matrix_entry($E4,[3,2],10); + + Context()->texStrings; + BEGIN_TEXT + The original matrix and its row reduced echelon form: + \[ $A \sim $Arref. \] + $BR + The original matrix with rows switched, multiplied, or added together: + \[ $Aswitch, $Amult, $Aadd. \] + $BR + Some elementary matrices. + \[$E1, $E2, $E3, $E4\] + END_TEXT + Context()->normalStrings; + + COMMENT('MathObject version.'); + ENDDOCUMENT(); =back diff --git a/macros/PGmatrixmacros.pl b/macros/PGmatrixmacros.pl index 070646c63a..85698a4fce 100644 --- a/macros/PGmatrixmacros.pl +++ b/macros/PGmatrixmacros.pl @@ -4,7 +4,7 @@ =head1 NAME - Matrix macros for the PG language + PGmatrixmacros.pl =head1 SYNPOSIS @@ -12,11 +12,22 @@ =head1 SYNPOSIS =head1 DESCRIPTION -Almost all of the macros in the file are very rough at best. The most useful is display_matrix. -Many of the other macros work with vectors and matrices stored as anonymous arrays. +Matrix macros for the PG language -Frequently it may be -more useful to use the Matrix objects defined RealMatrix.pm and Matrix.pm and the constructs listed there. +These macros are fairly old. The most useful is display_matrix and +its variants. + +Frequently it will be +most useful to use the MathObjects Matrix (defined in Value::Matrix.pm) +and Vector types which +have more capabilities and more error checking than the subroutines in +this file. These macros have no object orientation and +work with vectors and matrices +stored as perl anonymous arrays. + +There are also Matrix objects defined in +RealMatrix.pm and Matrix.pm but in almost all cases the +MathObjects Matrix types are preferable. =cut @@ -28,132 +39,57 @@ BEGIN sub _PGmatrixmacros_init { } -# this subroutine zero_check is not very well designed below -- if it is used much it should receive -# more work -- particularly for checking relative tolerance. More work needs to be done if this is -# actually used. - -sub zero_check{ - my $array = shift; - my %options = @_; - my $num = @$array; - my $i; - my $max = 0; my $mm; - for ($i=0; $i< $num; $i++) { - $mm = $array->[$i] ; - $max = abs($mm) if abs($mm) > $max; - } - my $tol = $options{tol}; - $tol = 0.01*$options{reltol}*$options{avg} if defined($options{reltol}) and defined $options{avg}; - $tol = .000001 unless defined($tol); - ($max <$tol) ? 1: 0; # 1 if the array is close to zero; -} -sub vec_dot{ - my $vec1 = shift; - my $vec2 = shift; - warn "vectors must have the same length" unless @$vec1 == @$vec2; # the vectors must have the same length. - my @vec1=@$vec1; - my @vec2=@$vec2; - my $sum = 0; - - while(@vec1) { - $sum += shift(@vec1)*shift(@vec2); - } - $sum; -} -sub proj_vec { - my $vec = shift; - warn "First input must be a column matrix" unless ref($vec) eq 'Matrix' and ${$vec->dim()}[1] == 1; - my $matrix = shift; # the matrix represents a set of vectors spanning the linear space - # onto which we want to project the vector. - warn "Second input must be a matrix" unless ref($matrix) eq 'Matrix' and ${$matrix->dim()}[1] == ${$vec->dim()}[0]; - $matrix * transpose($matrix) * $vec; -} - -sub vec_cmp{ #check to see that the submitted vector is a non-zero multiple of the correct vector - my $correct_vector = shift; - my %options = @_; - my $ans_eval = sub { - my $in = shift @_; - - my $ans_hash = new AnswerHash; - my @in = split("\0",$in); - my @correct_vector=@$correct_vector; - $ans_hash->{student_ans} = "( " . join(", ", @in ) . " )"; - $ans_hash->{correct_ans} = "( " . join(", ", @correct_vector ) . " )"; - - return($ans_hash) unless @$correct_vector == @in; # make sure the vectors are the same dimension - - my $correct_length = vec_dot($correct_vector,$correct_vector); - my $in_length = vec_dot(\@in,\@in); - return($ans_hash) if $in_length == 0; - - if (defined($correct_length) and $correct_length != 0) { - my $constant = vec_dot($correct_vector,\@in)/$correct_length; - my @difference = (); - for(my $i=0; $i < @correct_vector; $i++ ) { - $difference[$i]=$constant*$correct_vector[$i] - $in[$i]; - } - $ans_hash->{score} = zero_check(\@difference); - - } else { - $ans_hash->{score} = 1 if vec_dot(\@in,\@in) == 0; - } - $ans_hash; - - }; - - $ans_eval; -} ############ =head4 display_matrix - Usage \{ display_matrix( [ [1, '\(\sin x\)'], [ans_rule(5), 6] ]) \} - \{ display_matrix($A, align=>'crvl') \} - \[ \{ display_matrix_mm($A) \} \] - \[ \{ display_matrix_mm([ [1, 3], [4, 6] ]) \} \] - - display_matrix produces a matrix for display purposes. It checks whether - it is producing LaTeX output, or if it is displaying on a web page in one - of the various modes. The input can either be of type Matrix, Value::Matrix (mathobject) - or a reference to an array. - - Entries can be numbers, Fraction objects, bits of math mode, or answer - boxes. An entire row can be replaced by the string 'hline' to produce - a horizontal line in the matrix. - - display_matrix_mm functions similarly, except that it should be inside - math mode. display_matrix_mm cannot contain answer boxes in its entries. - Entries to display_matrix_mm should assume that they are already in - math mode. - - Both functions take an optional alignment string, similar to ones in - LaTeX tabulars and arrays. Here c for centered columns, l for left - flushed columns, and r for right flushed columns. - - The alignment string can also specify vertical rules to be placed in the - matrix. Here s or | denote a solid line, d is a dashed line, and v - requests the default vertical line. This can be set on a system-wide - or course-wide basis via the variable $defaultDisplayMatrixStyle, and - it can default to solid, dashed, or no vertical line (n for none). - - The matrix has left and right delimiters also specified by - $defaultDisplayMatrixStyle. They can be parentheses, square brackets, - braces, vertical bars, or none. The default can be overridden in - an individual problem with optional arguments such as left=>"|", or - right=>"]". - - You can specify an optional argument of 'top_labels'=> ['a', 'b', 'c']. - These are placed above the columns of the matrix (as is typical for - linear programming tableau, for example). The entries will be typeset - in math mode. - - Top labels require a bit of care. For image modes, they look better - with display_matrix_mm where it is all one big image, but they work with - display_matrix. With tth, you pretty much have to use display_matrix - since tth can't handle the TeX tricks used to get the column headers - up there if it gets the whole matrix at once. + Usage + \{ display_matrix( [ [1, '\(\sin x\)'], [ans_rule(5), 6] ]) \} + \{ display_matrix($A, align=>'crvl') \} + \[ \{ display_matrix_mm($A) \} \] + \[ \{ display_matrix_mm([ [1, 3], [4, 6] ]) \} \] + +display_matrix produces a matrix for display purposes. It checks whether +it is producing LaTeX output, or if it is displaying on a web page in one +of the various modes. The input can either be of type Matrix, Value::Matrix (mathobject) +or a reference to an array. + +Entries can be numbers, Fraction objects, bits of math mode, or answer +boxes. An entire row can be replaced by the string 'hline' to produce +a horizontal line in the matrix. + +display_matrix_mm functions similarly, except that it should be inside +math mode. display_matrix_mm cannot contain answer boxes in its entries. +Entries to display_matrix_mm should assume that they are already in +math mode. + +Both functions take an optional alignment string, similar to ones in +LaTeX tabulars and arrays. Here c for centered columns, l for left +flushed columns, and r for right flushed columns. + +The alignment string can also specify vertical rules to be placed in the +matrix. Here s or | denote a solid line, d is a dashed line, and v +requests the default vertical line. This can be set on a system-wide +or course-wide basis via the variable $defaultDisplayMatrixStyle, and +it can default to solid, dashed, or no vertical line (n for none). + +The matrix has left and right delimiters also specified by +$defaultDisplayMatrixStyle. They can be parentheses, square brackets, +braces, vertical bars, or none. The default can be overridden in +an individual problem with optional arguments such as left=>"|", or +right=>"]". + +You can specify an optional argument of 'top_labels'=> ['a', 'b', 'c']. +These are placed above the columns of the matrix (as is typical for +linear programming tableau, for example). The entries will be typeset +in math mode. + +Top labels require a bit of care. For image modes, they look better +with display_matrix_mm where it is all one big image, but they work with +display_matrix. With tth, you pretty much have to use display_matrix +since tth can't handle the TeX tricks used to get the column headers +up there if it gets the whole matrix at once. =cut @@ -692,6 +628,7 @@ sub mbox { =head4 ra_flatten_matrix Usage: ra_flatten_matrix($A) + returns: [a11, a12,a21,a22] where $A is a matrix object The output is a reference to an array. The matrix is placed in the array by iterating @@ -715,9 +652,17 @@ sub ra_flatten_matrix{ \@array; } -# This subroutine is probably obsolete and not generally useful. It was patterned after the APL -# constructs for multiplying matrices. It might come in handy for non-standard multiplication of -# of matrices (e.g. mod 2) for indice matrices. + +=head4 apl_matrix_mult() + + # This subroutine is probably obsolete and not generally useful. + # It was patterned after the APL + # constructs for multiplying matrices. It might come in handy + # for non-standard multiplication of + # of matrices (e.g. mod 2) for indice matrices. + +=cut + sub apl_matrix_mult{ my $ra_a= shift; my $ra_b= shift; @@ -763,9 +708,11 @@ sub make_matrix{ =head4 create2d_matrix -This can be a useful method for quickly entering small matrices by hand. --MEG +This can be a useful method for quickly entering small matrices by hand. + --MEG - create2d_matrix("1 2 4, 5 6 8"); + create2d_matrix("1 2 4, 5 6 8"); or + create2d_matrix("1 2 4; 5 6 8"); produces the anonymous array [[1,2,4],[5,6,8] ] @@ -775,11 +722,42 @@ =head4 create2d_matrix sub create2d_matrix { my $string = shift; - my @rows = split("\\s*,\\s*",$string); + my @rows = split("\\s*[,;]\\s*",$string); @rows = map {[split("\\s", $_ )]} @rows; [@rows]; } + +=head2 convert_to_array_ref { + + $output_matrix = convert_to_array_ref($input_matrix) + +Converts a MathObject matrix (ref($input_matrix eq 'Value::Matrix') +or a MatrixReal1 matrix (ref($input_matrix eq 'Matrix')to +a reference to an array (e.g [[4,6],[3,2]]). +This adaptor allows all of the LinearProgramming.pl subroutines to be used with +MathObject arrays. + +$mathobject_matrix->value outputs an array (usually an array of array references) so placing it inside +square bracket produces and array reference (of array references) which is what lp_display_mm() is +seeking. + +=cut + +sub convert_to_array_ref { + my $input = shift; + if (ref($input) eq 'Value::Matrix' ) { + $input = [$input->value]; + } elsif (ref($input) eq 'Matrix' ) { + $input = $input->array_ref; + } elsif (ref($input) =~/ARRAY/) { + # no change to input value + } else { + WARN_MESSAGE("This does not appear to be a matrix "); + } + $input; +} + =head4 check_matrix_from_ans_box_cmp An answer checker factory built on create2d_matrix. This still needs @@ -796,7 +774,6 @@ sub check_matrix_from_ans_box_cmp{ my $string_matrix_cmp = sub { $string = shift @_; my $studentMatrix; - # eval { $studentMatrix = Matrix(create2d_matrix($string)); die "I give up";}; #caught by op_mask $studentMatrix = Matrix(create2d_matrix($string)); die "I give up"; # main::DEBUG_MESSAGE(ref($studentMatrix). "$studentMatrix with error "); # errors are returned as warnings. Can't seem to trap them. @@ -815,67 +792,100 @@ sub check_matrix_from_ans_box_cmp{ } -=head2 convert_to_array_ref { - $output_matrix = convert_to_array_ref($input_matrix) +=head4 zero_check (deprecated -- use MathObjects matrices and vectors) -Converts a MathObject matrix (ref($input_matrix eq 'Value::Matrix') -or a MatrixReal1 matrix (ref($input_matrix eq 'Matrix')to -a reference to an array (e.g [[4,6],[3,2]]). -This adaptor allows all of the Linear Programming subroutines to be used with -MathObject arrays. + # this subroutine zero_check is not very well designed below -- if it is used much it should receive + # more work -- particularly for checking relative tolerance. More work needs to be done if this is + # actually used. -$mathobject_matrix->value outputs an array (usually an array of array references) so placing it inside -square bracket produces and array reference (of array references) which is what lp_display_mm() is -seeking. +=cut + +sub zero_check{ + my $array = shift; + my %options = @_; + my $num = @$array; + my $i; + my $max = 0; my $mm; + for ($i=0; $i< $num; $i++) { + $mm = $array->[$i] ; + $max = abs($mm) if abs($mm) > $max; + } + my $tol = $options{tol}; + $tol = 0.01*$options{reltol}*$options{avg} if defined($options{reltol}) and defined $options{avg}; + $tol = .000001 unless defined($tol); + ($max <$tol) ? 1: 0; # 1 if the array is close to zero; +} + +=head4 vec_dot() (deprecated -- use MathObjects vectors and matrices) + +sub vec_dot{ + my $vec1 = shift; + my $vec2 = shift; + warn "vectors must have the same length" unless @$vec1 == @$vec2; # the vectors must have the same length. + my @vec1=@$vec1; + my @vec2=@$vec2; + my $sum = 0; + + while(@vec1) { + $sum += shift(@vec1)*shift(@vec2); + } + $sum; +} + +=head4 proj_vect (deprecated -- use MathObjects vectors and matrices) =cut -sub convert_to_array_ref { - my $input = shift; - if (ref($input) eq 'Value::Matrix' ) { - $input = [$input->value]; - } elsif (ref($input) eq 'Matrix' ) { - $input = $input->array_ref; - } elsif (ref($input) =~/ARRAY/) { - # no change to input value - } else { - WARN_MESSAGE("This does not appear to be a matrix "); - } - $input; +sub proj_vec { + my $vec = shift; + warn "First input must be a column matrix" unless ref($vec) eq 'Matrix' and ${$vec->dim()}[1] == 1; + my $matrix = shift; # the matrix represents a set of vectors spanning the linear space + # onto which we want to project the vector. + warn "Second input must be a matrix" unless ref($matrix) eq 'Matrix' and ${$matrix->dim()}[1] == ${$vec->dim()}[0]; + $matrix * transpose($matrix) * $vec; +} + +=head4 vec_cmp (deprecated -- use MathObjects vectors and matrices) + +=cut + + +sub vec_cmp{ #check to see that the submitted vector is a non-zero multiple of the correct vector + my $correct_vector = shift; + my %options = @_; + my $ans_eval = sub { + my $in = shift @_; + + my $ans_hash = new AnswerHash; + my @in = split("\0",$in); + my @correct_vector=@$correct_vector; + $ans_hash->{student_ans} = "( " . join(", ", @in ) . " )"; + $ans_hash->{correct_ans} = "( " . join(", ", @correct_vector ) . " )"; + + return($ans_hash) unless @$correct_vector == @in; # make sure the vectors are the same dimension + + my $correct_length = vec_dot($correct_vector,$correct_vector); + my $in_length = vec_dot(\@in,\@in); + return($ans_hash) if $in_length == 0; + + if (defined($correct_length) and $correct_length != 0) { + my $constant = vec_dot($correct_vector,\@in)/$correct_length; + my @difference = (); + for(my $i=0; $i < @correct_vector; $i++ ) { + $difference[$i]=$constant*$correct_vector[$i] - $in[$i]; + } + $ans_hash->{score} = zero_check(\@difference); + + } else { + $ans_hash->{score} = 1 if vec_dot(\@in,\@in) == 0; + } + $ans_hash; + + }; + + $ans_eval; } -# sub format_answer{ -# my $ra_eigenvalues = shift; -# my $ra_eigenvectors = shift; -# my $functionName = shift; -# my @eigenvalues=@$ra_eigenvalues; -# my $size= @eigenvalues; -# my $ra_eigen = make_matrix( sub {my ($i,$j) = @_; ($i==$j) ? "e^{$eigenvalues[$j] t}": 0 }, $size,$size); -# my $out = qq! -# $functionName(t) =! . -# displayMatrix(apl_matrix_mult($ra_eigenvectors,$ra_eigen, -# 'times'=>sub{($_[0] and $_[1]) ? "$_[0]$_[1]" : ''}, -# 'plus'=>sub{ my $out = join("",@_); ($out) ?$out : '0' } -# ) ) ; -# $out; -# } -# sub format_vector_answer{ -# my $ra_eigenvalues = shift; -# my $ra_eigenvectors = shift; -# my $functionName = shift; -# my @eigenvalues=@$ra_eigenvalues; -# my $size= @eigenvalues; -# my $ra_eigen = make_matrix( sub {my ($i,$j) = @_; ($i==$j) ? "e^{$eigenvalues[$j] t}": 0 }, $size,$size); -# my $out = qq! -# $functionName(t) =! . -# displayMatrix($ra_eigenvectors)."e^{$eigenvalues[0] t}" ; -# $out; -# } -# sub format_question{ -# my $ra_matrix = shift; -# my $out = qq! y'(t) = ! . displayMatrix($B). q! y(t)! -# -# } 1; diff --git a/macros/PGmorematrixmacros.pl b/macros/PGmorematrixmacros.pl index b2af761786..8ed22573a9 100644 --- a/macros/PGmorematrixmacros.pl +++ b/macros/PGmorematrixmacros.pl @@ -5,9 +5,24 @@ BEGIN # set the prefix used for arrays. our $ArRaY = $main::PG->{ARRAY_PREFIX}; +=head2 NAME + + macros/PGmorematrixmacros.pl + +=cut + + sub _PGmorematrixmacros_init{} -sub random_inv_matrix { ## Builds and returns a random invertible \$row by \$col matrix. +=head4 random_inv_matrix + +## Builds and returns a random invertible \$row by \$col matrix. + +=cut + + +sub random_inv_matrix { +## Builds and returns a random invertible \$row by \$col matrix. warn "Usage: \$new_matrix = random_inv_matrix(\$rows,\$cols)" if (@_ != 2); @@ -54,16 +69,29 @@ sub random_diag_matrix{ ## Builds and returns a random diagonal \$n by \$n matri return $D; } +=head4 swap_rows ($matrix, $row1, $row2) + + (deprecated use MathObject Matrix instead) + +$matrix is assumed to be a RealMatrix1 object. +It is better to use MathObject Matrices and row swap mechanisms +from MatrixReduce.pl instead. + +=cut + + sub swap_rows{ warn "Usage: \$new_matrix = swap_rows(\$matrix,\$row1,\$row2);" if (@_ != 3); my $matrix = $_[0]; my ($i,$j) = ($_[1],$_[2]); + warn "Error: Rows to be swapped must exist!" if ($i>@$matrix or $j >@$matrix); warn "Warning: Swapping the same row is pointless" - if ($i==$j); + if ($i==$j); + my $cols = @{$matrix->[0]}; my $B = new Matrix(@$matrix,$cols); foreach my $k (1..$cols){ @@ -73,6 +101,16 @@ sub swap_rows{ return $B; } +=head4 row_mult ($matrix, $scaler, $row) + + (deprecated use MathObject Matrix instead) + +$matrix is assumed to be a RealMatrix1 object. +It is better to use MathObject Matrices and row swap mechanisms +from MatrixReduce.pl instead. + +=cut + sub row_mult{ warn "Usage: \$new_matrix = row_mult(\$matrix,\$scalar,\$row);" @@ -88,6 +126,16 @@ sub row_mult{ return $B; } +sub linear_combo($matrix, $scalar, $row1, $row2) + + (deprecated use MathObject Matrix instead) + +Adds a multiple of row1 to row2. + +$matrix is assumed to be a RealMatrix1 object. +It is better to use MathObject Matrices and subroutines +from MatrixReduce.pl instead. + sub linear_combo{ warn "Usage: \$new_matrix = linear_combo(\$matrix,\$scalar,\$row1,\$row2);" @@ -106,6 +154,15 @@ sub linear_combo{ return $B; } + +=head2 + +These should be compared to similar subroutines made later in +MatrixCheckers.pl + + +=cut + =head3 basis_cmp() Compares a list of vectors by finding the change of coordinate matrix @@ -378,6 +435,8 @@ sub compare_basis { =head2 vec_list_string +(this is mostly obsolete. One should use MathObject Vectors instead. ) + This is a check_syntax type method (in fact I borrowed some of that method's code) for vector input. The student needs to enter vectors like: [1,0,0],[1,2,3],[0,9/sqrt(10),1/sqrt(10)] Each entry can contain functions and operations and the usual math constants (pi and e). @@ -503,8 +562,14 @@ sub vec_list_string{ $rh_ans; } + + + =head5 ans_array_filter + (this filter is not necessary when using MathObjects. It may someday be useful + again if the AnswerEvaluator pipeline is used to its fullest extent. ) + This filter was created to get, format, and evaluate each entry of the ans_array and ans_array_extension answer entry methods. Running this filter is necessary to get all the entries out of the answer hash. Each entry is evaluated and the resulting number is put in the display for student answer @@ -616,6 +681,20 @@ sub ans_array_filter{ } +=head3 + +The following subroutines, meant to be used with MatrixReal1 type matrices, are +deprecated. In general you should use the MathObject Matrix type and the +checking methods in MatrixCheckers.pl + + are_orthogonal_vecs($vec_ref, %opts) + is_diagonal($matrix, %opts) + are_unit_vecs($vec_ref, %opts) + display_correct_vecs($vec_ref, %opts) + vec_solution_cmp($vec,%opts) + filter: compare_vec_solution($rh_ans,%opts); + +=cut sub are_orthogonal_vecs{ my ($vec_ref , %opts) = @_; diff --git a/macros/tableau.pl b/macros/tableau.pl index 90e2e16cf4..b571ff3bc5 100755 --- a/macros/tableau.pl +++ b/macros/tableau.pl @@ -6,180 +6,248 @@ ##### From gage_matrix_ops # 2014_HKUST_demo/templates/setSequentialWordProblem/bill_and_steve.pg:"gage_matrix_ops.pl", -=head3 Matrix extraction mechanisms +=head1 NAME - matrix_column_slice (was matrix_from_matrix_cols) + macros/tableau.pl - matrix_row_slice (was matrix_from_matrix_rows) + + +=head2 DESCRIPTION + + # We're going to have several types + # MathObject Matrices Value::Matrix + # tableaus form John Jones macros + # MathObject tableaus + # Containing an matrix $A coefficients for constraint + # A vertical vector $b for constants for constraints + # A horizontal vector $c for coefficients for objective function + # A vertical vector corresponding to the value $z of the objective function + # dimensions $n problem vectors, $m constraints = $m slack variables + # A basis Value::Set -- positions for columns which are independent and + # whose associated variables can be determined + # uniquely from the parameter variables. + # The non-basis (parameter) variables are set to zero. + # + # state variables (assuming parameter variables are zero or when given parameter variables) + # create the methods for updating the various containers + # create the method for printing the tableau with all its decorations + # possibly with switches to turn the decorations on and off. + + +The structure of the tableau is: + + ----------------------------------------------------------- + | | | | | + | A | S | 0 | b | + | | | | | + ---------------------------------------------- + | -c | 0 | 1 | 0 | + ---------------------------------------------- + Matrix A, the constraint matrix is n=num_problem_vars by m=num_slack_vars + Matrix S, the slack variables is m by m + Matrix b, the constraint constants is n by 1 + Matrix c, the objective function coefficients matrix is 1 by n or 2 by n + The next to the last column holds z or objective value + z(...x^i...) = c_i* x^i (Einstein summation convention) + FIXME: allow c to be a 2 by n matrix so that you can do phase1 calculations easily + + +=cut - matrix_extract_submatrix (was matrix_from_submatrix) + +=head2 Package main + +=cut + + +=item tableauEquivalence + + ANS( $tableau->cmp(checker=>tableauEquivalence()) ); - matrix_extract_rows + # Note: it is important to include the () at the end of tableauEquivalence + + # tableauEquivalence is meant to compare two matrices up to + # reshuffling the rows and multiplying each row by a constant. + # E.g. equivalent up to multiplying on the left by a permuation matrix + # or a (non-uniformly constant) diagonal matrix + +=cut + +=item get_tableau_variable_values + + Parameters: ($MathObjectMatrix_tableau, $MathObjectSet_basis) + Returns: ARRAY or ARRAY_ref - matrix_extract_columns +Returns the solution variables to the tableau assuming +that the parameter (non-basis) variables +have been set to zero. It returns a list in +array context and a reference to +an array in scalar context. + +=item lp_basis_pivot - matrix_columns_to_List + Parameters: ($old_tableau,$old_basis,$pivot) + Returns: ($new_tableau, Set($new_basis),\@statevars) - matrix_rows_to_List +=item linebreak_at_commas + + Parameters: () + Return: -Many of these duplicate methods of Value::Matrix -- refactor. + Useage: + $foochecker = $constraints->cmp()->withPostFilter( + linebreak_at_commas() + ); + +Replaces commas with line breaks in the latex presentations of the answer checker. +Used most often when $constraints is a LinearInequality math object. + + + +=head2 Package tableau + +=item Tableau->new(A=>Matrix, b=>Vector or Matrix, c=>Vector or Matrix) + + A => undef, # original constraint matrix MathObjectMatrix + b => undef, # constraint constants Vector or MathObjectMatrix 1 by n + c => undef, # coefficients for objective function Vector or MathObjectMatrix 1 by n + obj_row => undef, # contains the negative of the coefficients of the objective function. + z => undef, # value for objective function + n => undef, # dimension of problem variables (columns in A) + m => undef, # dimension of slack variables (rows in A) + S => undef, # square m by m matrix for slack variables + M => undef, # matrix of consisting of all original columns and all + rows except for the objective function row + obj_col_num => undef, + basis => undef, # list describing the current basis columns corresponding + to determined variables. + current_basis_matrix => undef, # square invertible matrix + corresponding to the current basis columns + + current_constraint_matrix=>undef, # the current version of [A | S] + current_b, # the current version of the constraint constants b + # current_basis_matrix # (should be new name for B above + # # a square invertible matrix corresponding to the + # # current basis columns) + # flag indicating the column (1 or n+m+1) for the objective value + constraint_labels => undef, (not sure if this remains relevant after pivots) + problem_var_labels => undef, + slack_var_labels => undef, + +=item $self->current_tableau + Parameters: () + Returns: A MathObjectMatrix_tableau + +This represents the current version of the tableau + +=item $self->objective_row + Parameters: () + Returns: + +=item $self->basis + Parameter: ARRAY or ARRAY_ref or () + Returns: MathObject_list + + FiXME -- this should accept a MathObject_List (or MO_Set?) + +=head3 Package Tableau (eventually package Matrix?) + +=item $self->row_slice + + Parameter: @slice or \@slice + Return: MathObject matrix + +=item $self->extract_rows + + Parameter: @slice or \@slice + Return: two dimensional array ref + +=item extract_rows_to_list + + Parameter: @slice or \@slice + Return: MathObject List of row references + +=item $self->extract_columns + + Parameter: @slice or \@slice + Return: two dimensional array ref + +=item $self->column_slice + + Parameter: @slice or \@slice + Return: MathObject Matrix + +=item $self->extract_columns_to_list + + Parameter: @slice or \@slice + Return: MathObject List of Matrix references ? + +=item $self->submatrix + + Parameter:(rows=>\@row_slice,columns=>\@column_slice) + Return: MathObject matrix + =cut +=head3 References: + +MathObject Matrix methods: L +MathObject Contexts: L +CPAN RealMatrix docs: L + +More references: L + +=cut + sub _tableau_init {}; # don't reload this file -package main; -sub matrix_column_slice{ - matrix_from_matrix_cols(@_); -} -sub matrix_from_matrix_cols { - my $M = shift; # a MathObject matrix_columns - my($n,$m) = $M->dimensions; - my @slice = @_; - if (ref($slice[0]) =~ /ARRAY/) { # handle array reference - @slice = @{$slice[0]}; - } - @slice = @slice?@slice : (1..$m); - my @columns = map {$M->column($_)->transpose->value} @slice; - #create the chosen columns as rows - # then transform to array_refs. - Matrix(@columns)->transpose; #transpose and return an n by m matrix (2 dim) -} -sub matrix_row_slice{ - matrix_from_matrix_rows(@_); -} +# loadMacros("tableau_main_subroutines.pl"); -sub matrix_from_matrix_rows { - my $M = shift; # a MathObject matrix_columns - my($n,$m) = $M->dimensions; - my @slice = @_; - if (ref($slice[0]) =~ /ARRAY/) { # handle array reference - @slice = @{$slice[0]}; - } - @slice = @slice? @slice : (1..$n); # the default is the whole matrix. - # DEBUG_MESSAGE("row slice in matrix from rows is @slice"); - my @rows = map {[$M->row($_)->value]} @slice; - #create the chosen columns as rows - # then transform to array_refs. - Matrix([@rows]); # insure that it is still an n by m matrix (2 dim) -} -sub matrix_extract_submatrix { - matrix_from_submatrix(@_); -} -sub matrix_from_submatrix { - my $M=shift; - return undef unless ref($M) =~ /Value::Matrix/; - my %options = @_; - my($n,$m) = $M->dimensions; - my $row_slice = ($options{rows})?$options{rows}:[1..$m]; - my $col_slice = ($options{columns})?$options{columns}:[1..$n]; - # DEBUG_MESSAGE("ROW SLICE", join(" ", @$row_slice)); - # DEBUG_MESSAGE("COL SLICE", join(" ", @$col_slice)); - my $M1 = matrix_from_matrix_rows($M,@$row_slice); - # DEBUG_MESSAGE("M1 - matrix from rows) $M1"); - return matrix_from_matrix_cols($M1, @$col_slice); -} -sub matrix_extract_rows { - my $M =shift; - my @slice = @_; - if (ref($slice[0]) =~ /ARRAY/) { # handle array reference - @slice = @{$slice[0]}; - } elsif (@slice == 0) { # export all rows to List - @slice = ( 1..(($M->dimensions)[0]) ); - } - return map {$M->row($_)} @slice ; -} +=head4 Subroutines added to the main:: Package -sub matrix_rows_to_list { - List(matrix_extract_rows(@_)); -} -sub matrix_columns_to_list { - List(matrix_extract_columns(@_) ); -} -sub matrix_extract_columns { - my $M =shift; # Add error checking - my @slice = @_; - if (ref($slice[0]) =~ /ARRAY/) { # handle array reference - @slice = @{$slice[0]}; - } elsif (@slice == 0) { # export all columns to an array - @slice = 1..($M->dimensions->[1]); + +=cut + +=item tableauEquivalence + + + +=cut + +sub tableauEquivalence { + return sub { + my ($correct, $student, $ansHash) = @_; + # DEBUG_MESSAGE("executing tableau equivalence"); + # convert matrices to arrays of row references + my @rows1 = matrix_extract_rows($correct); + my @rows2 = matrix_extract_rows($student); + # compare the rows as lists with each row being compared as + # a parallel Vector (i.e. up to multiples) + my $score = List(@rows1)->cmp( checker => + sub { + my ($listcorrect,$liststudent,$listansHash,$nth,$value)=@_; + my $listscore = Vector($listcorrect)->cmp(parallel=>1) + ->evaluate(Vector($liststudent))->{score}; + return $listscore; + } + )->evaluate(List(@rows2))->{score}; + return $score; } - return map {$M->column($_)} @slice; -} + } +=item linebreak_at_commas + + Useage: + + $foochecker = $constraints->cmp()->withPostFilter( + linebreak_at_commas() + ); +=cut -######################## -############## -# get_tableau_variable_values -# -# Calculates the values of the basis variables of the tableau, -# assuming the parameter variables are 0. -# -# Usage: get_tableau_variable_values($MathObjectMatrix_tableau, $MathObjectSet_basis) -# -# feature request -- for tableau object -- allow specification of non-zero parameter variables -sub get_tableau_variable_values { - my $mat = shift; # a MathObject matrix - my $basis =shift; # a MathObject set - # FIXME - # type check ref($mat)='Matrix'; ref($basis)='Set'; - # or check that $mat has dimensions, element methods; and $basis has a contains method - my ($n, $m) = $mat->dimensions; - @var = (); - #DEBUG_MESSAGE( "start new matrix"); - foreach my $j (1..$m-2) { # the last two columns of the tableau are object variable and constants - if (not $basis->contains($j)) { - #DEBUG_MESSAGE( "j= $j not in basis"); - $var[$j-1]=0; next; # non-basis variables (parameters) are set to 0. - - } else { - foreach my $i (1..$n-1) { # the last row is the objective function - # if this is a basis column there should be only one non-zero element(the pivot) - if ( not $mat->element($i,$j) == 0 ) { # should this have ->value????? - $var[$j-1] = ($mat->element($i,$m)/$mat->element($i,$j))->value; - #DEBUG_MESSAGE("i=$i j=$j var = $var[$j-1] "); - next; - } - - } - } - } # element($n, $m-1) is the coefficient of the objective value. - # this last variable is the value of the objective function - push @var , ($mat->element($n,$m)/$mat->element($n,$m-1))->value; - - @var; -} -#### Test -- assume matrix is this -# 1 2 1 0 0 | 0 | 3 -# 4 5 0 1 0 | 0 | 6 -# 7 8 0 0 1 | 0 | 9 -# -1 -2 0 0 0 | 1 | 10 # objective row -# and basis is {3,4,5} (start columns with 1) -# $n= 4; $m = 7 -# $x1=0; $x2=0; $x3=s1=3; $x4=s2=6; $x5=s3=9; w=10=objective value -# -# - -#################################### -# -# Cover for lp_pivot which allows us to use a set object for the new and old basis - -sub lp_basis_pivot { - my ($old_tableau,$old_basis,$pivot) = @_; # $pivot is a Value::Point - my $new_tableau= lp_clone($old_tableau); - main::lp_pivot($new_tableau, $pivot->extract(1)-1,$pivot->extract(2)-1); - my $new_matrix = Matrix($new_tableau); - my ($n,$m) = $new_matrix->dimensions; - my $param_size = $m-$n -1; #n=constraints+1, #m = $param_size + $constraints +2 - my $new_basis = ( $old_basis - ($pivot->extract(1)+$param_size) + ($pivot->extract(2)) )->sort; - my @statevars = get_tableau_variable_values($new_matrix, $new_basis); - return ( $new_tableau, Set($new_basis),\@statevars); #FIXME -- force to set (from type Union) to insure that ->data is an array and not a string. -} - - sub linebreak_at_commas { return sub { my $ans=shift; @@ -194,37 +262,40 @@ sub linebreak_at_commas { $ans; }; } -# Useage -# $foochecker = $constraints->cmp()->withPostFilter( -# linebreak_at_commas() -# ); - - -# We're going to have several types -# MathObject Matrices Value::Matrix -# tableaus form John Jones macros -# MathObject tableaus -# Containing an matrix $A coefficients for constraint -# A vertical vector $b for constants for constraints -# A horizontal vector $c for coefficients for objective function -# A vertical vector $P for the value of the objective function -# dimensions $n problem vectors, $m constraints = $m slack variables -# A basis Value::Set -- positions for columns which are independent and -# whose associated variables can be determined -# uniquely from the parameter variables. -# The non-basis (parameter) variables are set to zero. -# -# state variables (assuming parameter variables are zero or when given parameter variables) -# create the methods for updating the various containers -# create the method for printing the tableau with all its decorations -# possibly with switches to turn the decorations on and off. - - -### End gage_matrix_ops include + +=item linebreak_at_commas + + Useage: + + lop_display($tableau, align=>'cccc|cc|c|c', toplevel=>[qw(x1,x2,x3,x4,s1,s2,P,b)]) + +Pretty prints the output of a matrix as a LOP with separating labels and +variable labels. + +=cut +sub lop_display { + my $tableau = shift; + %options = @_; + $options{alignment} = ($options{alignment})? $options{alignment}:"|ccccc|cc|c|c|"; + @toplevel = (); + if (exists( ($options{toplevel})) ) { + @toplevel = @{$options{toplevel}}; + $toplevel[0]=[$toplevel[0],headerrow=>1, midrule=>1]; + } + @matrix = $tableau1->current_tableau->value; + $last_row = $#matrix; # last row is objective coefficients + $matrix[$last_row-1]->[0]=[$matrix[$last_row-1]->[0],midrule=>1]; + $matrix[$last_row]->[0]=[$matrix[$last_row]->[0],midrule=>1]; + DataTable([[@toplevel],@matrix],align=>$options{alignment}); +} + + ################################################## package Tableau; our @ISA = qw(Value::Matrix Value); +sub class {"Matrix"}; + sub _Matrix { # can we just import this? Value::Matrix->new(@_); } @@ -235,7 +306,7 @@ sub new { my $tableau = { A => undef, # constraint matrix MathObjectMatrix b => undef, # constraint constants Vector or MathObjectMatrix 1 by n - c => undef, # coefficients for objective function Vector or MathObjectMatrix 1 by n + c => undef, # coefficients for objective function MathObjectMatrix 1 by n or 2 by n matrix obj_row => undef, # contains the negative of the coefficients of the objective function. z => undef, # value for objective function n => undef, # dimension of problem variables (columns in A) @@ -244,10 +315,14 @@ sub new { basis => undef, # list describing the current basis columns corresponding to determined variables. B => undef, # square invertible matrix corresponding to the current basis columns M => undef, # matrix of consisting of all columns and all rows except for the objective function row - obj_col_num => undef, # flag indicating the column (1 or n+m+1) for the objective value + current_constraint_matrix=>undef, + current_coeff=>undef, + current_b => undef, + obj_col_index => undef, # an array reference indicating the columns (e.g 1 or n+m+1) for the objective value or values constraint_labels => undef, problem_var_labels => undef, slack_var_labels => undef, + @_, }; bless $tableau, $class; @@ -255,25 +330,35 @@ sub new { return $tableau; } + +# the following are used to construct the tableau +# initialize +# assemble_matrix +# objective_row + sub initialize { $self= shift; unless (ref($self->{A}) =~ /Value::Matrix/ && ref($self->{b}) =~ /Value::Vector|Value::Matrix/ && ref($self->{c}) =~ /Value::Vector|Value::Matrix/){ - main::WARN_MESSAGE("Error: Required inputs: Tableau(A=> Matrix, b=>Vector, c=>Vector)"); - return; + Value::Error ("Error: Required inputs for creating tableau:\n". + "Tableau(A=> Matrix, b=>ColumnVector or Matrix, c=>Vector or Matrix)"); } my ($m, $n)=($self->{A}->dimensions); $self->{n}=$self->{n}//$n; $self->{m}=$self->{m}//$m; # main::DEBUG_MESSAGE("m $m, n $n"); - $self->{S} = Value::Matrix->I(4); + $self->{S} = Value::Matrix->I($m); $self->{basis} = [($n+1)...($n+$m)] unless ref($self->{basis})=~/ARRAY/; + my @rows = $self->assemble_matrix; - #main::DEBUG_MESSAGE("rows", @rows); - $self->{M} = _Matrix([@rows]); - $self->{B} = $self->{M}->submatrix(rows=>[1..($self->{m})],columns=>$self->{basis}); - $self->{obj_row} = _Matrix($self->objective_row()); + # main::DEBUG_MESSAGE("rows", map {ref($_)?$_->value :$_} map {@$_} @rows); + $self->{M} = _Matrix([@rows]); #original matrix + $self->{current_constraint_matrix}= $self->{M}; + $self->{data}= $self->{M}->data; + $self->{current_basis_matrix} = $self->{M}->submatrix(rows=>[1..($self->{m})],columns=>$self->{basis}); + $self->{obj_row} = _Matrix(@{$self->objective_row()}); + $self->basis($self->basis()->value); return(); } @@ -282,79 +367,458 @@ sub assemble_matrix { my @rows =(); my $m = $self->{m}; my $n = $self->{n}; + # sanity check for b; + if (ref($self->{b}) =~/Vector/) { + # replace by n by 1 matrix + $self->{b}=Value::Matrix->new([[$self->{b}->value]])->transpose; + } + my ($constraint_rows, $constraint_cols) = $self->{b}->dimensions; + unless ($constraint_rows== $m and $constraint_cols == 1 ) { + Value::Error("constraint matrix b is $constraint_rows by $constraint_cols but should + be $m by 1 to match the constraint matrix A "); + } + foreach my $i (1..$m) { my @current_row=(); foreach my $j (1..$n) { - push @current_row, $self->{A}->element($i, $j); + push @current_row, $self->{A}->element($i, $j)->value; } foreach my $j (1..$m) { - push @current_row, $self->{S}->element($i,$j); # slack variables + push @current_row, $self->{S}->element($i,$j)->value; # slack variables } - push @current_row, 0, $self->{b}->data->[$i-1]; # obj column and constant column + push @current_row, 0, $self->{b}->row($i)->value; # obj column and constant column push @rows, [@current_row]; } return @rows; # these are the matrices A | S | obj | b - # the final row describing the objective function is not in this + # the final row describing the objective function + # is not in this part of the matrix } sub objective_row { my $self = shift; + # sanity check for objective row + + my @last_row=(); - push @last_row, ( -($self->{c}) )->value; - foreach my $i (1..($self->{m})) { push @last_row, 0 }; - push @last_row, 1, 0; + push @last_row, ( -($self->{c}) )->value; # add the negative coefficients of the obj function + foreach my $i (1..($self->{m})) { push @last_row, 0 }; # add 0s for the slack variables + push @last_row, 1, 0; # add the 1 for the objective value and 0 for the initial value return \@last_row; } +=item current_tableau + + Useage: + $MathObjectmatrix = $self->current_tableau + $MathObjectmatrix = $self->current_tableau(3,4) #updates basis to (3,4) + +Returns the current constraint matrix, including the constraint constants AND the +row containing the objective function coefficients. + +=cut +sub current_tableau { + my $self = shift; + my @basis = @_; + if (@basis) { + $self->basis(@basis); + } + # find the coefficients associated with the basis columns + my $c_B = $self->{obj_row}->extract_columns($self->{basis} ); + my $c_B2 = Value::Vector->new([ map {$_->value} @$c_B]); +# my $current_constraint_matrix = $self->{current_constraint_matrix}; +# my $B = $self->{current_basis_matrix}; + #main::DEBUG_MESSAGE( "basis: current_constraint_matrix $current_constraint_matrix"); + #main::DEBUG_MESSAGE( "current_tableau: matrix is $current_constraint_matrix"); + #main::DEBUG_MESSAGE( "current_basis matrix: $B"); +# my $correction_coeff = ($c_B2*($self->{current_constraint_matrix}) )->row(1); + # subtract the correction coefficients from the obj_row + # this essentially extends Gauss reduction applied to the obj_row +# my $obj_row_normalized = abs($self->{current_basis_matrix}->det->value)*$self->{obj_row}; + #main::DEBUG_MESSAGE(" normalized obj row ",$obj_row_normalized->value); + #main::DEBUG_MESSAGE(" correction coeff ", $correction_coeff->value); +# my $current_coeff = $obj_row_normalized-$correction_coeff ; +# $self->{current_coeff}= $current_coeff; + + #main::DEBUG_MESSAGE("subtract these two ", (($self->{current_basis_matrix}->det) *$self->{obj_row}), " | ", ($c_B*$current_tableau)->dimensions); + #main::DEBUG_MESSAGE("all coefficients", join('|', $self->{obj_row}->value ) ); + #main::DEBUG_MESSAGE("current coefficients", join('|', @current_coeff) ); + #main::DEBUG_MESSAGE("type of $self->{basis}", ref($self->{basis}) ); + #main::DEBUG_MESSAGE("current basis",join("|", @{$self->{basis}})); + #main::DEBUG_MESSAGE("CURRENT STATE ", $current_tableau); + return _Matrix( @{$self->{current_constraint_matrix}->extract_rows},$self->{current_coeff} ); + #return( $self->{current_coeff} ); +} + sub basis { my $self = shift; #update basis + # basis is stored as an ARRAY reference. + # basis is exported as a list + # FIXME should basis be sorted? + Value::Error( "call basis as a method") unless ref($self)=~/Tableau/; my @input = @_; return Value::List->new($self->{basis}) unless @input; #return basis if no input my $new_basis; if (ref( $input[0]) =~/ARRAY/) { $new_basis=$input[0]; - } else { + } elsif (ref( $input[0]) =~/List|Set/){ + $new_basis = $input[0]->value; + } else { # input is assumed to be an array $new_basis = \@input; } - $self->{basis}= $new_basis; - $self->{B} = $self->{M}->submatrix(rows=>[1..($self->{m})],columns=>$self->{basis}); - return Value::List->new($self->{basis}); -} - -sub current_state { - my $Badj = ($self->{B}->det) * ($self->{B}->inverse); - my $current_tableau = $Badj * $self->{M}; # the A | S | obj | b - $self->{current_tableau}=$current_tableau; - # find the coefficients associated with the basis columns + $self->{basis}= $new_basis; # this should always be an ARRAY + WARN_MESSAGE("basis $new_basis was not stored as an array reference") + unless ref($new_basis)=~/ARRAY/; + + + $self->{current_basis_matrix} = $self->{M}->submatrix(rows=>[1..($self->{m})],columns=>$self->{basis}); + my $B = $self->{current_basis_matrix}; + #main::DEBUG_MESSAGE("basis: B is $B" ); + my $Badj = abs($self->{current_basis_matrix}->det->value) * ($self->{current_basis_matrix}->inverse); + my $M = $self->{M}; + my ($row_dim, $col_dim) = $M->dimensions; + my $current_constraint_matrix = $Badj*$M; my $c_B = $self->{obj_row}->extract_columns($self->{basis} ); my $c_B2 = Value::Vector->new([ map {$_->value} @$c_B]); - my $correction_coeff = ($c_B2*$current_tableau )->row(1); - # subtract the correction coefficients from the obj_row - # this is essentially extends Gauss reduction applied to the obj_row - my $obj_row_normalized = ($self->{B}->det) *$self->{obj_row}; + my $correction_coeff = ($c_B2*($current_constraint_matrix) )->row(1); + my $obj_row_normalized = abs($self->{current_basis_matrix}->det->value)*$self->{obj_row}; my $current_coeff = $obj_row_normalized-$correction_coeff ; + # updates + $self->{current_basis_matrix}= $B; + $self->{data} = $current_constraint_matrix->data; + $self->{current_constraint_matrix} = $current_constraint_matrix; $self->{current_coeff}= $current_coeff; + $self->{current_b} = $current_constraint_matrix->column($col_dim); + # the A | S | obj | b + # main::DEBUG_MESSAGE( "basis: current_constraint_matrix $current_constraint_matrix ". + # ref($self->{current_constraint_matrix}) ); + # main::DEBUG_MESSAGE("basis self ",ref($self), "---", ref($self->{basis})); + return Value::List->new($self->{basis}); +} - #main::DEBUG_MESSAGE("subtract these two ", (($self->{B}->det) *$self->{obj_row}), " | ", ($c_B*$current_tableau)->dimensions); - #main::DEBUG_MESSAGE("all coefficients", join('|', $self->{obj_row}->value ) ); - #main::DEBUG_MESSAGE("current coefficients", join('|', @current_coeff) ); - #main::DEBUG_MESSAGE("type of $self->{basis}", ref($self->{basis}) ); - #main::DEBUG_MESSAGE("current basis",join("|", @{$self->{basis}})); - #main::DEBUG_MESSAGE("CURRENT STATE ", $current_state); - return _Matrix( @{$current_tableau->extract_rows},$self->{current_coeff} ); - #return( $self->{current_coeff} ); + +=item find_next_basis + + ($row, $col,$optimum,$unbounded) = $self->find_next_basis (max/min, obj_row_number) + +In phase 2 of the simplex method calculates the next basis. +$optimum or $unbounded is set +if the process has found on the optimum value, or the column +$col gives a certificate of unboundedness. + + +=cut + + +sub find_next_basis { + my $self = shift; + my $max_or_min = shift; + my $obj_row_number = shift; + my ( $row_index, $col_index, $optimum, $unbounded)= + $self->find_next_pivot($max_or_min, $obj_row_number); + my $flag; + my $basis; + if ($optimum or $unbounded) { + $basis=$self->basis(); + } else { + $flag = ''; + $basis =$self->find_next_basis_from_pivot($row_index,$col_index); + + } + return( $basis->value, $optimum,$unbounded ); + +} + +=item find_next_pivot + + ($row, $col,$optimum,$unbounded) = $self->find_next_pivot (max/minm obj_row_number) + +This is used in phase2 so the possible outcomes are only $optimum and $unbounded. +$infeasible is not possible. Use the lowest index strategy to find the next pivot +point. This calls find_pivot_row and find_pivot_column. $row and $col are undefined if +either $optimum or $unbounded is set. + +=cut + +sub find_next_pivot { + my $self = shift; + my $max_or_min = shift; + my $obj_row_number =shift; + + # sanity check max or min in find pivot column + my ($col_index, $value, $row_index, $optimum, $unbounded) = ('','','',''); + ($col_index, $value, $optimum) = $self->find_pivot_column($max_or_min, $obj_row_number); +# main::DEBUG_MESSAGE("find_next_pivot: col: $col_index, value: $value opt: $optimum "); + return ( $row_index, $col_index, $optimum, $unbounded) if $optimum; + ($row_index, $value, $unbounded) = $self->find_pivot_row($col_index); +# main::DEBUG_MESSAGE("find_next pivot row: $row_index, value: $value unbound: $unbounded"); + return($row_index, $col_index, $optimum, $unbounded); +} + + + +=item find_next_basis_from_pivot + + List(basis) = $self->find_next_basis (row_index, col_index) + +Calculate the next basis from the current basis given the pivot position. + +=cut + +sub find_next_basis_from_pivot { + my $self = shift; + my $row_index = shift; + my $col_index =shift; + # sanity check max or min in find pivot column + my $basis = main::Set($self->{basis}); + my ($leaving_col_index, $value) = $self->find_leaving_column($row_index); + $basis = main::Set( $basis - main::Set($leaving_col_index) + main::Set($col_index)); + # main::DEBUG_MESSAGE( "basis is $basis, leaving index $leaving_col_index + # entering index is $col_index"); + #$basis = [$basis->value, main::Real($col_index)]; + return ($basis); +} + + + +=item find_pivot_column + + ($index, $value, $optimum) = $self->find_pivot_column (max/min, obj_row_number) + +This finds the left most obj function coefficient that is negative (for maximizing) +or positive (for minimizing) and returns the value and the index. Only the +index is really needed for this method. The row number is included because there might +be more than one objective function in the table (for example when using +the Auxiliary method in phase1 of the simplex method.) If there is no coefficient +of the appropriate sign then the $optimum flag is set and $index and $value +are undefined. + +=cut + +sub find_pivot_column { + my $self = shift; + my $max_or_min = shift; + my $obj_row_index = shift; + # sanity check + unless ($max_or_min =~ /max|min/) { + Value::Error( "The optimization method must be + 'max' or 'min'. |$max_or_min| is not defined."); + } + my $obj_row_matrix = $self->{current_coeff}; + #FIXME $obj_row_matrix is this a 1 by n or an n dimensional matrix?? + my ($obj_col_dim) = $obj_row_matrix->dimensions; + my $obj_row_dim = 1; + $obj_col_dim=$obj_col_dim-2; + #sanity check row + if (not defined($obj_row_index) ) { + $obj_row_index = 1; + } elsif ($obj_row_index<1 or $obj_row_index >$obj_row_dim){ + Value::Error( "The choice for the objective row $obj_row_index is out of range."); + } + #FIXME -- make sure objective row is always a two dimensional matrix, often with one row. + + + my @obj_row = @{$obj_row_matrix->extract_rows($obj_row_index)}; + my $index = -1; + my $optimum = 1; + my $value = undef; +# main::DEBUG_MESSAGE(" coldim: $obj_col_dim , row: $obj_row_index obj_matrix: $obj_row_matrix ".ref($obj_row_matrix) ); +# main::DEBUG_MESSAGE(" \@obj_row ", join(' ', @obj_row ) ); + for (my $k=1; $k<=$obj_col_dim; $k++) { +# main::DEBUG_MESSAGE("find pivot column: k $k, " .$obj_row_matrix->element($k)->value); + if ( ($obj_row_matrix->element($k) < 0 and $max_or_min eq 'max') or + ($obj_row_matrix->element($k) > 0 and $max_or_min eq 'min') ) { + $index = $k; #memorize index + $value = $obj_row_matrix->element($k); + $optimum = 0; + last; # found first coefficient with correct sign + } + } + return ($index, $value, $optimum); +} + +=item find_pivot_row + + ($index, $value, $unbounded) = $self->find_pivot_row(col_number) + +Compares the ratio $b[$j]/a[$j, $col_number] and chooses the smallest +non-negative entry. It assumes that we are in phase2 of simplex methods +so that $b[j]>0; If all entries are negative (or infinity) then +the $unbounded flag is set and returned and the $index and $value +quantities are undefined. + +=cut + +sub find_pivot_row { + my $self = shift; + my $column_index = shift; + my ($row_dim, $col_dim) = $self->{M}->dimensions; + $col_dim = $col_dim-2; # omit the obj_value and constraint columns + # sanity check column_index + unless (1<=$column_index and $column_index <= $col_dim) { + Value::Error( "Column index must be between 1 and $col_dim" ); + } + # main::DEBUG_MESSAGE("dim = ($row_dim, $col_dim)"); + my $value = undef; + my $index = -1; + my $unbounded = 1; + for (my $k=1; $k<=$row_dim; $k++) { + my $m = $self->{current_constraint_matrix}->element($k,$column_index); + # main::DEBUG_MESSAGE(" m[$k,$column_index] is ", $m->value); + next if $m <=0; + my $b = $self->{current_b}->element($k,1); + # main::DEBUG_MESSAGE(" b[$k] is ", $b->value); + # main::DEBUG_MESSAGE("finding pivot row in column $column_index, row: $k ", ($b/$m)->value); + if ( not defined($value) or $b/$m < $value) { + $value = $b/$m; + $index = $k; # memorize index + $unbounded = 0; + } + } + return( $index, $value, $unbounded); +} + + + + +=item find_leaving_column + + ($index, $value) = $self->find_leaving_column(obj_row_number) + +Finds the non-basis column with a non-zero entry in the given row. When +called with the pivot row number this index gives the column which will +be removed from the basis while the pivot col number gives the basis +column which will become a parameter column. + +=cut + +sub find_leaving_column { + my $self = shift; + my $row_index = shift; + my ($row_dim,$col_dim) = $self->{current_constraint_matrix}->dimensions; + $col_dim= $col_dim - 1; # both problem and slack variables are included + # but not the constraint column or the obj_value column(s) (the latter are zero) + # FIXME + #sanity check row index; + unless (1<=$row_index and $row_index <= $row_dim) { + Value::Error("The row number must be between 1 and $row_dim" ); + } + my $basis = main::Set($self->{basis}); + my $index = 0; + my $value = undef; + foreach my $k (1..$col_dim) { + next unless $basis->contains(main::Set($k)); + $m_ik = $self->{current_constraint_matrix}->element($row_index, $k); + next unless $m_ik !=0; + $index = $k; # memorize index + $value = $m_ik; + last; + } + return( $index, $value); } +=item find_short_cut_row + + ($index, $value, $feasible)=$self->find_short_cut_row + +Find the most negative entry in the constraint column vector $b. If all entries +are positive then the tableau represents a feasible point, $feasible is set to 1 +and $index and $value are undefined. + +=cut + +sub find_short_cut_row { + my $self = shift; +# my @b = map {$_->value} @{ $self->{b}->column(1)->data }; #get raw numbers +# main::DEBUG_MESSAGE("b is ", join(" ", @b )); + my ($row_dim, $col_dim) = $self->{b}->dimensions; +# main::DEBUG_MESSAGE($self->{b}->dimensions); +# main::DEBUG_MESSAGE("dimensions of b: $row_dim, $col_dim "); + my $col_index = 1; + my $index = undef; + my $value = undef; + my $feasible = 1; + for (my $k=1; $k<=$row_dim; $k++) { + my $b_k1 = $self->{current_b}->element($k,$col_index); + # main::DEBUG_MESSAGE("b[$k] = $b_k1"); + next if $b_k1>=0; #skip positive entries; + $index =$k; + $value = $b_k1; + $feasible = 0; + last; + } + return ( $index, $value, $feasible); +} + +=item find_short_cut_column + + ($index, $value, $infeasible) = $self->find_short_cut_column(row number) + +Find the left most negative entry in the specified row. If all coefficients are +positive then the tableau represents an infeasible LOP, the $infeasible flag is set, +and the $index and $value are undefined. + +=cut + +sub find_short_cut_column { + my $self = shift; + my $row_index = shift; + my ($row_dim,$col_dim) = $self->{M}->dimensions; + $col_dim = $col_dim - 1; # omit constraint column + # FIXME to adjust for additional obj_value columns + #sanity check row index + unless (1<= $row_index and $row_index <= $row_dim) { + Value::Error("The row must be between 1 and $row_dim"); + } + my $index = undef; + my $value = undef; + my $infeasible = 1; + for (my $k = 1; $k<=$col_dim; $k++ ) { + my $m_ik = $self->{current_constraint_matrix}->element($row_index, $k); + # main::DEBUG_MESSAGE( "in M: ($row_index, $k) contains $m_ik"); + next if $m_ik >=0; + $index = $k; + $value = $m_ik; + $infeasible = 0; + last; + } + return( $index, $value, $infeasible); +} + +=item find_next_short_cut_pivot + + ($row, $col, $feasible, $infeasible) = $self->find_next_short_cut_pivot + + +Following the short-cut algorithm this chooses the next pivot by choosing the row +with the most negative constraint constant entry (top most first in case of tie) and +then the left most negative entry in that row. + +The process stops with either $feasible=1 (state variables give a feasible point for the +constraints) or $infeasible=1 (a row in the tableau shows that the LOP has empty domain.) + +=cut +=item find_next_short_cut_basis + +=cut + +# eventually these routines should be included in the Value::Matrix +# module? package Value::Matrix; -sub _Matrix { +sub _Matrix { Value::Matrix->new(@_); } + +sub row_slice { + my $self = shift; + @slice = @_; + return _Matrix( $self->extract_rows(@slice) ); +} sub extract_rows { - $self = shift; + my $self = shift; my @slice = @_; if (ref($slice[0]) =~ /ARRAY/) { # handle array reference @slice = @{$slice[0]}; @@ -363,9 +827,12 @@ sub extract_rows { } return [map {$self->row($_)} @slice ]; #prefer to pass references when possible } - -sub extract_columns { - $self = shift; +sub column_slice { + my $self = shift; + return _Matrix( $self->extract_columns(@_) )->transpose; # matrix is built as rows then transposed. +} +sub extract_columns { + my $self = shift; my @slice = @_; if (ref($slice[0]) =~ /ARRAY/) { # handle array reference @slice = @{$slice[0]}; @@ -386,19 +853,6 @@ sub extract_columns_to_list { Value::List->new($self->extract_columns(@_) ); } -sub column_slice { - $self = shift; - return _Matrix( $self->extract_columns(@_) )->transpose; # matrix is built as rows then transposed. -} - - -sub row_slice { - $self = shift; - @slice = @_; - return _Matrix( $self->extract_rows(@slice) ); -} - - sub submatrix { my $self = shift; my %options = @_; diff --git a/macros/tableau_main_subroutines.pl b/macros/tableau_main_subroutines.pl new file mode 100644 index 0000000000..8eb3249854 --- /dev/null +++ b/macros/tableau_main_subroutines.pl @@ -0,0 +1,186 @@ +# subroutines included into the main:: package. + +package main; + +sub isMatrix { + my $m = shift; + return ref($m) =~/Value::Matrix/i; +} +sub matrix_column_slice{ + matrix_from_matrix_cols(@_); +} +sub matrix_from_matrix_cols { + my $M = shift; # a MathObject matrix_columns + my($n,$m) = $M->dimensions; + my @slice = @_; + if (ref($slice[0]) =~ /ARRAY/) { # handle array reference + @slice = @{$slice[0]}; + } + @slice = @slice?@slice : (1..$m); + my @columns = map {$M->column($_)->transpose->value} @slice; + #create the chosen columns as rows + # then transform to array_refs. + Matrix(@columns)->transpose; #transpose and return an n by m matrix (2 dim) +} +sub matrix_row_slice{ + matrix_from_matrix_rows(@_); +} + +sub matrix_from_matrix_rows { + my $M = shift; # a MathObject matrix_columns + unless (isMatrix($M)){ + WARN_MESSAGE( "matrix_from_matrix_rows: |".ref($M)."| or |$M| is not a MathObject Matrix"); + return undef; + } + + my($n,$m) = $M->dimensions; + my @slice = @_; + if (ref($slice[0]) =~ /ARRAY/) { # handle array reference + @slice = @{$slice[0]}; + } + @slice = @slice? @slice : (1..$n); # the default is the whole matrix. + # DEBUG_MESSAGE("row slice in matrix from rows is @slice"); + my @rows = map {[$M->row($_)->value]} @slice; + #create the chosen columns as rows + # then transform to array_refs. + Matrix([@rows]); # insure that it is still an n by m matrix (2 dim) +} + +sub matrix_extract_submatrix { + matrix_from_submatrix(@_); +} +sub matrix_from_submatrix { + my $M=shift; + unless (isMatrix($M)){ + warn( "matrix_from_submatrix: |".ref($M)."| or |$M| is not a MathObject Matrix"); + return undef; + } + + my %options = @_; + my($n,$m) = $M->dimensions; + my $row_slice = ($options{rows})?$options{rows}:[1..$m]; + my $col_slice = ($options{columns})?$options{columns}:[1..$n]; + # DEBUG_MESSAGE("ROW SLICE", join(" ", @$row_slice)); + # DEBUG_MESSAGE("COL SLICE", join(" ", @$col_slice)); + my $M1 = matrix_from_matrix_rows($M,@$row_slice); + # DEBUG_MESSAGE("M1 - matrix from rows) $M1"); + return matrix_from_matrix_cols($M1, @$col_slice); +} +sub matrix_extract_rows { + my $M =shift; + unless (isMatrix($M)){ + WARN_MESSAGE( "matrix_extract_rows: |".ref($M)."| or |$M| is not a MathObject Matrix"); + return undef; + } + + my @slice = @_; + if (ref($slice[0]) =~ /ARRAY/) { # handle array reference + @slice = @{$slice[0]}; + } elsif (@slice == 0) { # export all rows to List + @slice = ( 1..(($M->dimensions)[0]) ); + } + return map {$M->row($_)} @slice ; +} + +sub matrix_rows_to_list { + List(matrix_extract_rows(@_)); +} +sub matrix_columns_to_list { + List(matrix_extract_columns(@_) ); +} +sub matrix_extract_columns { + my $M =shift; # Add error checking + unless (isMatrix($M)){ + WARN_MESSAGE( "matrix_extract_columns: |".ref($M)."| or |$M| is not a MathObject Matrix"); + return undef; + } + + my @slice = @_; + if (ref($slice[0]) =~ /ARRAY/) { # handle array reference + @slice = @{$slice[0]}; + } elsif (@slice == 0) { # export all columns to an array + @slice = 1..($M->dimensions->[1]); + } + return map {$M->column($_)} @slice; +} + + + +######################## +############## +# get_tableau_variable_values +# +# Calculates the values of the basis variables of the tableau, +# assuming the parameter variables are 0. +# +# Usage: ARRAY = get_tableau_variable_values($MathObjectMatrix_tableau, $MathObjectSet_basis) +# +# feature request -- for tableau object -- allow specification of non-zero parameter variables +sub get_tableau_variable_values { + my $mat = shift; # a MathObject matrix + unless (isMatrix($mat)){ + WARN_MESSAGE( "get_tableau_variable_values: |".ref($mat)."| or |$mat| is not a MathObject Matrix"); + return Matrix([0]); + } + my $basis =shift; # a MathObject set + # FIXME + # type check ref($mat)='Matrix'; ref($basis)='Set'; + # or check that $mat has dimensions, element methods; and $basis has a contains method + my ($n, $m) = $mat->dimensions; + @var = (); + #DEBUG_MESSAGE( "start new matrix"); + foreach my $j (1..$m-2) { # the last two columns of the tableau are object variable and constants + if (not $basis->contains($j)) { + # DEBUG_MESSAGE( "j= $j not in basis"); # set the parameter values to zero + $var[$j-1]=0; next; # non-basis variables (parameters) are set to 0. + + } else { + foreach my $i (1..$n-1) { # the last row is the objective function + # if this is a basis column there should be only one non-zero element(the pivot) + if ( $mat->element($i,$j)->value != 0 ) { # should this have ->value????? + $var[$j-1] = ($mat->element($i,$m)/($mat->element($i,$j))->value); + # DEBUG_MESSAGE("i=$i j=$j var = $var[$j-1] "); # calculate basis variable value + next; + } + + } + } + } # element($n, $m-1) is the coefficient of the objective value. + # this last variable is the value of the objective function + # check for division by zero + if ($mat->element($n,$m-1)->value != 0 ) { + push @var , ($mat->element($n,$m)/$mat->element($n,$m-1))->value; + } else { + push @var , '.666'; + } + return wantarray ? @var : \@var; # return either array or reference to an array +} +#### Test -- assume matrix is this +# 1 2 1 0 0 | 0 | 3 +# 4 5 0 1 0 | 0 | 6 +# 7 8 0 0 1 | 0 | 9 +# -1 -2 0 0 0 | 1 | 10 # objective row +# and basis is {3,4,5} (start columns with 1) +# $n= 4; $m = 7 +# $x1=0; $x2=0; $x3=s1=3; $x4=s2=6; $x5=s3=9; w=10=objective value +# +# + +#################################### +# +# Cover for lp_pivot which allows us to use a set object for the new and old basis + +sub lp_basis_pivot { + my ($old_tableau,$old_basis,$pivot) = @_; # $pivot is a Value::Point + my $new_tableau= lp_clone($old_tableau); + # lp_pivot has 0 based indices + main::lp_pivot($new_tableau, $pivot->extract(1)-1,$pivot->extract(2)-1); + # lp_pivot pivots in place + my $new_matrix = Matrix($new_tableau); + my ($n,$m) = $new_matrix->dimensions; + my $param_size = $m-$n -1; #n=constraints+1, #m = $param_size + $constraints +2 + my $new_basis = ( $old_basis - ($pivot->extract(1)+$param_size) + ($pivot->extract(2)) )->sort; + my @statevars = get_tableau_variable_values($new_matrix, $new_basis); + return ( $new_tableau, Set($new_basis),\@statevars); + # force to set (from type Union) to insure that ->data is an array and not a string. +} diff --git a/t/matrix_tableau_tests/matrix_test3.pg b/t/matrix_tableau_tests/matrix_test3.pg new file mode 100644 index 0000000000..4c8892b662 --- /dev/null +++ b/t/matrix_tableau_tests/matrix_test3.pg @@ -0,0 +1,74 @@ +############################################## +DOCUMENT(); + +loadMacros( + "PGstandard.pl", # Standard macros for PG language + "MathObjects.pl", + "parserLinearInequality.pl", + "PGML.pl", + "quickMatrixEntry.pl", + "tableau.pl", + "PGmatrixmacros.pl", + "LinearProgramming.pl", + #"source.pl", # allows code to be displayed on certain sites. + "PGcourse.pl", +); + +############################################## + +Context("Matrix"); # need Matrix context to allow string input into Matrix. +$zeroLevelTol = 1E-03; +$zeroLevel = 1E-03; +$tolType = 'relative', # "absolute"; #relative +$tolerance = .001; + +Context()->flags->set( + tolType => $tolType, + tolerance =>$tolerance, + zeroLevel=>$zeroLevel, + zeroLevelTol=>$zeroLevelTol +); +our $context=Context(); + + +BEGIN_PGML +Test to see if the tolerance flags are being respected. +zeroLevelTol = [$zeroLevelTol] +zeroLevel = [$zeroLevel] +tolType = [$tolType] +tolerance = [$tolerance] +END_PGML +#Construct a small test matrix. +$m = Matrix("[[3,6,7],[2,0,8],[4,6,21],[-6,4,0]]"); +$m2 = Matrix("[[3,6,7],[1.9999,.0001,8],[4,6,21],[-6,4,1E-5]]"); + +BEGIN_PGML +m is [$m] +and m2 is [$m2] + +Enter [$m2] +[@MATRIX_ENTRY_BUTTON($m2)@]* +[@ ANS($m2->cmp), $m2->ans_array(15) @]* + + +Are these tolerance leveals being ignored? +m and m2 are equal = [@ $m == $m2 @] + +END_PGML + +BEGIN_TEXT +This is in the BEGIN_TEXT/END_TEXT section: $BR + +m is $m $BR +and m2 is $m2 $BR + +Enter $m +\{ MATRIX_ENTRY_BUTTON($m) \} +\{ ANS($m->cmp), $m->ans_array \} +$PAR +FIXME -- these tolerance leveals are being ignored +m and m2 are equal = \{ $m == $m2 \} + +END_TEXT + +ENDDOCUMENT(); \ No newline at end of file diff --git a/t/matrix_tableau_tests/matrix_test4.pg b/t/matrix_tableau_tests/matrix_test4.pg new file mode 100644 index 0000000000..82109a548f --- /dev/null +++ b/t/matrix_tableau_tests/matrix_test4.pg @@ -0,0 +1,78 @@ +############################################## +DOCUMENT(); + +loadMacros( + "PGstandard.pl", # Standard macros for PG language + "MathObjects.pl", + "parserLinearInequality.pl", + "PGML.pl", + "quickMatrixEntry.pl", + "tableau.pl", + "PGmatrixmacros.pl", + "LinearProgramming.pl", + #"source.pl", # allows code to be displayed on certain sites. + "PGcourse.pl", +); + +############################################## + +Context("Matrix"); # need Matrix context to allow string input into Matrix. +$zeroLevelTol = 1E-03; +$zeroLevel = 1E-03; +$tolType = 'relative'; # "absolute"; #relative +$tolerance = .001; + +Context("Matrix"); + +Context()->flags->set( + tolType => $tolType, + tolerance =>$tolerance, + zeroLevel=>$zeroLevel, + zeroLevelTol=>$zeroLevelTol +); + + +BEGIN_PGML +Test to see if the tolerance flags are being respected +when inserted into answer checkers +zeroLevelTol = [$zeroLevelTol] +zeroLevel = [$zeroLevel] +tolType = [$tolType] +tolerance = [$tolerance] +END_PGML + +#Construct a small test matrix. + +$m = Matrix("[[3,6,7],[2,0,8],[4,6,21],[-6,4,0]]"); +$m2 = Matrix("[[3,6,7],[1.9999,.0001,8],[4,6,21],[-6,4,1E-5]]"); + +BEGIN_PGML +m is [$m] +and m2 is [$m2] + +Enter [$m2] +[@MATRIX_ENTRY_BUTTON($m2)@]* +[@ ANS($m2->cmp), $m2->ans_array(15) @]* + + +Are these tolerance leveals being ignored? +m and m2 are equal = [@ $m == $m2 @] + +END_PGML + +BEGIN_TEXT +This is in the BEGIN_TEXT/END_TEXT section: $BR + +m is $m $BR +and m2 is $m2 $BR + +Enter $m +\{ MATRIX_ENTRY_BUTTON($m) \} +\{ ANS($m->cmp), $m->ans_array \} +$PAR +Are these tolerance levels are being ignored? +m and m2 are equal = \{ $m == $m2 \} + +END_TEXT + +ENDDOCUMENT(); \ No newline at end of file diff --git a/t/matrix_tableau_tests/tableau-test6.pg b/t/matrix_tableau_tests/tableau-test6.pg new file mode 100644 index 0000000000..547a762075 --- /dev/null +++ b/t/matrix_tableau_tests/tableau-test6.pg @@ -0,0 +1,177 @@ +DOCUMENT(); +loadMacros( +"PGstandard.pl", +"MathObjects.pl", +"parserPopUp.pl", +"unionLists.pl", +"MatrixReduce.pl", +#"AppletObjects.pl", +"PGessaymacros.pl", +"PGmatrixmacros.pl", +"LinearProgramming.pl", +"parserLinearInequality.pl", +"quickMatrixEntry.pl", +#"scaffold.pl", +"tableau.pl", +#"gage_matrix_ops.pl", +"PGinfo.pl", +"source.pl", +"PGcourse.pl", +); + +TEXT(beginproblem()); +TEXT($BEGIN_ONE_COLUMN); +$showPartialCorrectAnswers = 1; + +INITIALIZE_QUICK_MATRIX_ENTRY(); + +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + + +#Hack to prevent domain conflict in answer. +Context()->variables->add(p1=>'Real',p2=>'Real'); +$objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); +# why can't the formula be defined within context "linearInequality"? + +Context("LinearInequality"); +Context()->variables->add(p1=>'Real',p2=>'Real'); +Context()->strings->add("Essay Answer" =>{}); +Context()->strings->add('Minimize'=>{},'Maximize'=>{}, "?"=>{}); +Context()->strings->add('Yes'=>{},'No'=>{}); +Context()->flags->set({tolType=>"absolute",tolerance=>.001}); +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +# $objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); + +our $context=Context(); + + +sub update_tableau { + my $original_matrix = shift; + my $basis = shift; + # $basis is a Set object + # last row is numbered $basis_size+1; + my $basis_size=@{$basis->data}; + my $basis_matrix = $original_matrix->submatrix(rows=>[1..$basis_size],columns=>($basis->data)); + my $context=Context(); + my $det = 1; #Real($basis_matrix->det); + $det = Real($basis_matrix->det); + $det = $det->value; + $det=Real(130000); #Real(1300000); + DEBUG_MESSAGE("det = $det ", ref($det)); + my $basis2= Matrix($basis_matrix); + #DEBUG_MESSAGE(($basis2->det)*($basis2->inverse) *$original_matrix->row_slice(1..$basis_size), $PAR); + my $obj_coeff = - ($original_matrix->row($basis_size+1)); + #obj function coeff corresponding to basis variables/columns + $cB = Matrix($obj_coeff->column_slice($basis->data)); + # get A | S |b portion of tableau and multiply by det(B)B^-1 + my $current_matrix= ($det)*($basis_matrix->inverse)*($original_matrix->row_slice(1..$basis_size)); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + return $current_tableau; +} +# the full tableau has one more row -- the objective function +sub display_tableau_state { + my ($original_matrix,$tableau, $matrix, $basis,$pivot)= @_; + $basis = $basis->sort; + my $basis_size = @{$basis->data} ; + my $basis_matrix = matrix_from_submatrix($original_matrix, rows=>[1..$basis_size],columns=>$basis->data); + my $normalized_tableau = $matrix; + my $reduced_matrix = matrix_from_matrix_rows($matrix,1..$basis_size); + warn "basis matrix is $basis_matrix, $reduced_matrix"; + my $normalized_reduced_matrix = ($basis_matrix->det)*$reduced_matrix; + # calculate current tableau by multiplying by the inverse of the basis + my $obj_coeff = - $original_matrix->row($basis_size+1); + my $cB = Matrix($obj_coeff->column_slice($basis->data)); + my $current_matrix = ($basis_matrix->det)*($basis_matrix->inverse)*matrix_from_matrix_rows($original_matrix,1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + $current_tableau = update_tableau($original_matrix, $basis); + + #DEBUG_MESSAGE("get state variables for $matrix using basis $basis"); + @statevars1 = get_tableau_variable_values($matrix, $basis); + # get z value + $statevars1 = ~~@statevars1; + $state = Matrix([[@statevars1]]); + + + return " + + pivot: \($pivot\) basis: \( $basis\) $PAR basis matrix \(". + display_matrix_mm($basis_matrix)."\) $PAR + tableau and normalized current tableau $PAR + \(" . lp_display_mm($tableau) . "\)\(".lp_display_mm($normalized_tableau). " \) $PAR + original matrix and current matrix calculated by Binverse*original $PAR + \(". + display_matrix_mm(matrix_from_matrix_rows($original_matrix,1..$basis_size)). + "\) \(". + display_matrix_mm($current_matrix). + "\) $PAR + objective coefficients = \( $obj_coeff \) $PAR + basis objective_coefficients = \($cB\) $PAR + new tableau \[". + lp_display_mm($current_tableau). + "\] $PAR + state: \($state)\) $PAR + " +} + + +#### problem starts here: + +Context($context); +$ra_matrix = [[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1,0,0,0,-$bill_profit], + [-1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, 0,1,0,0,-$steve_profit], + [1,0,0,0,0,0,0,0,1,0]]; + +# [0,$money_total,$time_total,1,1,0,0,1,0,0]]; +$original_matrix = Matrix($ra_matrix); + +$ra_matrix1 = lp_clone($ra_matrix); +Context()->texStrings; +BEGIN_TEXT +basis(6,7) $BR +\[\{lp_display_mm($ra_matrix1)\}\] $PAR +\{lp_pivot($ra_matrix1,0,0),''\}$BR pivot(1,1) basis(1,7)$PAR +\[\{lp_display_mm($ra_matrix1)\}\] $PAR +\{lp_pivot($ra_matrix1,1,1),''\}$BR pivot(2,2), basis(1,2) $PAR +\[\{lp_display_mm($ra_matrix1)\}\] $PAR +\{lp_pivot($ra_matrix1,0,2),''\}$BR pivot(1,3), basis(2,3) $PAR +\[\{lp_display_mm($ra_matrix1)\}\] $PAR +END_TEXT + + +$basis2_3_string = update_tableau($original_matrix, Set(2,3)); + + +BEGIN_TEXT +update tableau method $PAR +basis(6,7)$BR +\[ \{ update_tableau($original_matrix, Set(6,7)) \} \] $PAR +basis(1,7) $BR +\[ \{ update_tableau($original_matrix, Set(1,7)) \} \] $PAR +basis(1,2) $BR +\[ \{ update_tableau($original_matrix, Set(1,2)) \} \] $PAR +basis(2,4) $BR +\[ \{ update_tableau($original_matrix, Set(2,4) ) \} \] $PAR +basis(2,3) $BR +\[ $basis2_3_string \] $PAR +basis(2,3) $BR +\[ $basis2_3_string \] $PAR +\[ update_tableau($original_matrix, Set(2,3)\] +END_TEXT +Context()->normalStrings; +TEXT($END_ONE_COLUMN); +ENDDOCUMENT(); diff --git a/t/matrix_tableau_tests/tableau-test6a.pg b/t/matrix_tableau_tests/tableau-test6a.pg new file mode 100644 index 0000000000..d13e303684 --- /dev/null +++ b/t/matrix_tableau_tests/tableau-test6a.pg @@ -0,0 +1,186 @@ +DOCUMENT(); +loadMacros( +"PGstandard.pl", +"MathObjects.pl", +"parserPopUp.pl", +"unionLists.pl", +"MatrixReduce.pl", +#"AppletObjects.pl", +"PGessaymacros.pl", +"PGmatrixmacros.pl", +"LinearProgramming.pl", +"parserLinearInequality.pl", +"quickMatrixEntry.pl", +#"scaffold.pl", +"tableau.pl", +#"gage_matrix_ops.pl", +"PGinfo.pl", +"source.pl", +"PGcourse.pl", +); + +TEXT(beginproblem()); +TEXT($BEGIN_ONE_COLUMN); +$showPartialCorrectAnswers = 1; + +INITIALIZE_QUICK_MATRIX_ENTRY(); + +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + + +#Hack to prevent domain conflict in answer. +Context()->variables->add(p1=>'Real',p2=>'Real'); +$objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); +# why can't the formula be defined within context "linearInequality"? + +Context("LinearInequality"); +Context()->variables->add(p1=>'Real',p2=>'Real'); +Context()->strings->add("Essay Answer" =>{}); +Context()->strings->add('Minimize'=>{},'Maximize'=>{}, "?"=>{}); +Context()->strings->add('Yes'=>{},'No'=>{}); +Context()->flags->set({tolType=>"absolute",tolerance=>.001}); +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +# $objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); + +our $context=Context(); + +sub basis_size { + my $basis = shift; + return scalar( @{$basis->data} ); +} +sub basis_matrix { + my $original_matrix = shift; + my $basis = shift; + my $basis_size=basis_size($basis); + my $basis_matrix = $original_matrix->submatrix(rows=>[1..$basis_size],columns=>($basis->data)); + return $basis_matrix; +} + +sub update_tableau { + my $original_matrix = shift; + my $basis = shift; + # $basis is a Set object + # last row is numbered $basis_size+1; + my $basis_matrix = basis_matrix($original_matrix); + my $det = Real($basis_matrix->det); +# $det = Real($basis_matrix->det); +# $det = $det->value; +# $det=Real(130000); #Real(1300000); +# DEBUG_MESSAGE("det = $det ", ref($det)); +# my $basis2= Matrix($basis_matrix); +# #DEBUG_MESSAGE(($basis2->det)*($basis2->inverse) *$original_matrix->row_slice(1..$basis_size), $PAR); + my $obj_coeff = - ($original_matrix->row($basis_size+1)); + #obj function coeff corresponding to basis variables/columns + $cB = Matrix($obj_coeff->column_slice($basis->data)); + # get A | S |b portion of tableau and multiply by det(B)B^-1 + my $current_matrix= ($det)*($basis_matrix->inverse)*($original_matrix->row_slice(1..$basis_size)); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + return $current_tableau; +} +# the full tableau has one more row -- the objective function +sub display_tableau_state { + my ($original_matrix,$tableau, $matrix, $basis,$pivot)= @_; + $basis = $basis->sort; + my $basis_size = @{$basis->data} ; + my $basis_matrix = matrix_from_submatrix($original_matrix, rows=>[1..$basis_size],columns=>$basis->data); + my $normalized_tableau = $matrix; + my $reduced_matrix = matrix_from_matrix_rows($matrix,1..$basis_size); + warn "basis matrix is $basis_matrix, $reduced_matrix"; + my $normalized_reduced_matrix = ($basis_matrix->det)*$reduced_matrix; + # calculate current tableau by multiplying by the inverse of the basis + my $obj_coeff = - $original_matrix->row($basis_size+1); + my $cB = Matrix($obj_coeff->column_slice($basis->data)); + my $current_matrix = ($basis_matrix->det)*($basis_matrix->inverse)*matrix_from_matrix_rows($original_matrix,1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + $current_tableau = update_tableau($original_matrix, $basis); + + #DEBUG_MESSAGE("get state variables for $matrix using basis $basis"); + @statevars1 = get_tableau_variable_values($matrix, $basis); + # get z value + $statevars1 = ~~@statevars1; + $state = Matrix([[@statevars1]]); + + + return " + + pivot: \($pivot\) basis: \( $basis\) $PAR basis matrix \(". + display_matrix_mm($basis_matrix)."\) $PAR + tableau and normalized current tableau $PAR + \(" . lp_display_mm($tableau) . "\)\(".lp_display_mm($normalized_tableau). " \) $PAR + original matrix and current matrix calculated by Binverse*original $PAR + \(". + display_matrix_mm(matrix_from_matrix_rows($original_matrix,1..$basis_size)). + "\) \(". + display_matrix_mm($current_matrix). + "\) $PAR + objective coefficients = \( $obj_coeff \) $PAR + basis objective_coefficients = \($cB\) $PAR + new tableau \[". + lp_display_mm($current_tableau). + "\] $PAR + state: \($state)\) $PAR + " +} + + +#### problem starts here: + +Context($context); +$ra_matrix = [[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1,0,0,0,-$bill_profit], + [-1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, 0,1,0,0,-$steve_profit], + [1,0,0,0,0,0,0,0,1,0]]; + +# [0,$money_total,$time_total,1,1,0,0,1,0,0]]; +$original_matrix = Matrix($ra_matrix); + +$ra_matrix1 = lp_clone($ra_matrix); + +TEXT("no pivot det67 $det67$BR"); +$det67 = $det = basis_matrix($original_matrix,Set(6,7))->det->value; +$display1 = lp_display_mm($det67*Matrix($ra_matrix)); + +TEXT("det17 $det17$BR"); +$det17 = $det = basis_matrix($original_matrix,Set(1,7))->det; +$display2=lp_display_mm(Matrix( lp_pivot($ra_matrix1,0,0) )); + +TEXT("pivot (2,2), det12 $det12$BR"); +$det12 = $det = basis_matrix($original_matrix,Set(1,2))->det; +$display3 = lp_display_mm(Matrix( lp_pivot($ra_matrix1,1,1) )); + +TEXT("pivot (1,3) det23 $det23$BR"); +$det23 = $det = basis_matrix($original_matrix,Set(2,3))->det; +$display4 = lp_display_mm(Matrix( lp_pivot($ra_matrix1,0,2) )); + +Context()->texStrings; + +BEGIN_TEXT +tableau-test6a.pg $BR + +no pivot -- basis(6,7) $PAR +\[ $display1 \] $PAR + pivot(1,1) basis(1,7)$PAR +\[ $display2 \] $PAR +pivot(2,2), basis(1,2) $PAR +\[ $display3 \] $PAR +pivot(1,3), basis(2,3) $PAR +\[ $display4 \] $PAR +END_TEXT + +Context()->normalStrings; +TEXT($END_ONE_COLUMN); +ENDDOCUMENT(); diff --git a/t/matrix_tableau_tests/tableau-test6b.pg b/t/matrix_tableau_tests/tableau-test6b.pg new file mode 100644 index 0000000000..c0d95c1184 --- /dev/null +++ b/t/matrix_tableau_tests/tableau-test6b.pg @@ -0,0 +1,190 @@ +DOCUMENT(); +loadMacros( +"PGstandard.pl", +"MathObjects.pl", +"parserPopUp.pl", +"unionLists.pl", +"MatrixReduce.pl", +#"AppletObjects.pl", +"PGessaymacros.pl", +"PGmatrixmacros.pl", +"LinearProgramming.pl", +"parserLinearInequality.pl", +"quickMatrixEntry.pl", +#"scaffold.pl", +"tableau.pl", +#"gage_matrix_ops.pl", +"PGinfo.pl", +"source.pl", +"PGcourse.pl", +); + +TEXT(beginproblem()); +TEXT($BEGIN_ONE_COLUMN); +$showPartialCorrectAnswers = 1; + +INITIALIZE_QUICK_MATRIX_ENTRY(); + +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + + +#Hack to prevent domain conflict in answer. +Context()->variables->add(p1=>'Real',p2=>'Real'); +$objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); +# why can't the formula be defined within context "linearInequality"? + +Context("LinearInequality"); +Context()->variables->add(p1=>'Real',p2=>'Real'); +Context()->strings->add("Essay Answer" =>{}); +Context()->strings->add('Minimize'=>{},'Maximize'=>{}, "?"=>{}); +Context()->strings->add('Yes'=>{},'No'=>{}); +Context()->flags->set({tolType=>"absolute",tolerance=>.001}); +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +# $objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); + +our $context=Context(); + +sub basis_size { + my $basis = shift; + return scalar( @{$basis->data} ); +} +sub basis_matrix { + my $original_matrix = shift; + my $basis = shift; + my $basis_size=basis_size($basis); + my $basis_matrix = $original_matrix->submatrix(rows=>[1..$basis_size],columns=>($basis->data)); + return $basis_matrix; +} + +sub update_tableau { + my $original_matrix = shift; + my $basis = shift; #basis columns -- a Set + # $basis is a Set object + # last row is numbered $basis_size+1; + my $context= Context(); + #Context()->normalStrings; + my $basis_matrix = basis_matrix($original_matrix, $basis); + my $basis_size=basis_size($basis); + my $det = Real($basis_matrix->det)->value; # prevent tex output from det + my $obj_coeff = - ($original_matrix->row($basis_size+1)); + #obj function coeff corresponding to basis variables/columns + $cB = Matrix($obj_coeff->column_slice($basis->data)); + # get A | S |b portion of tableau and multiply by det(B)B^-1 + my $current_matrix= ($det)*($basis_matrix->inverse)*($original_matrix->row_slice(1..$basis_size)); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + #Context()->texStrings; # restore context + return $current_tableau; +} +# the full tableau has one more row -- the objective function +sub display_tableau_state { + my ($original_matrix,$tableau, $matrix, $basis,$pivot)= @_; + $basis = $basis->sort; + my $basis_size = @{$basis->data} ; + my $basis_matrix = matrix_from_submatrix($original_matrix, rows=>[1..$basis_size],columns=>$basis->data); + my $normalized_tableau = $matrix; + my $reduced_matrix = matrix_from_matrix_rows($matrix,1..$basis_size); + warn "basis matrix is $basis_matrix, $reduced_matrix"; + my $normalized_reduced_matrix = ($basis_matrix->det)*$reduced_matrix; + # calculate current tableau by multiplying by the inverse of the basis + my $obj_coeff = - $original_matrix->row($basis_size+1); + my $cB = Matrix($obj_coeff->column_slice($basis->data)); + my $current_matrix = ($basis_matrix->det)*($basis_matrix->inverse)*matrix_from_matrix_rows($original_matrix,1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + $current_tableau = update_tableau($original_matrix, $basis); + + #DEBUG_MESSAGE("get state variables for $matrix using basis $basis"); + @statevars1 = get_tableau_variable_values($matrix, $basis); + # get z value + $statevars1 = ~~@statevars1; + $state = Matrix([[@statevars1]]); + + + return " + + pivot: \($pivot\) basis: \( $basis\) $PAR basis matrix \(". + display_matrix_mm($basis_matrix)."\) $PAR + tableau and normalized current tableau $PAR + \(" . lp_display_mm($tableau) . "\)\(".lp_display_mm($normalized_tableau). " \) $PAR + original matrix and current matrix calculated by Binverse*original $PAR + \(". + display_matrix_mm(matrix_from_matrix_rows($original_matrix,1..$basis_size)). + "\) \(". + display_matrix_mm($current_matrix). + "\) $PAR + objective coefficients = \( $obj_coeff \) $PAR + basis objective_coefficients = \($cB\) $PAR + new tableau \[". + lp_display_mm($current_tableau). + "\] $PAR + state: \($state)\) $PAR + " +} + + +#### problem starts here: + +Context($context); +$ra_matrix = [[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1,0,0,0,-$bill_profit], + [-1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, 0,1,0,0,-$steve_profit], + [1,0,0,0,0,0,0,0,1,0]]; + +# [0,$money_total,$time_total,1,1,0,0,1,0,0]]; +$original_matrix = Matrix($ra_matrix); + +$ra_matrix1 = lp_clone($ra_matrix); + + +$string67 = lp_display_mm( update_tableau($original_matrix, Set(6,7)) ); +$string17 = lp_display_mm( update_tableau($original_matrix, Set(1,7)) ); +$string12 = lp_display_mm( update_tableau($original_matrix, Set(1,2)) ); +$string23 = lp_display_mm( update_tableau($original_matrix, Set(2,3)) ); +$string32 = lp_display_mm( update_tableau($original_matrix, Set(3,2)) ); + +Context()->texStrings; +BEGIN_TEXT +update tableau method tableau-test6b$PAR +basis(6,7)$BR +\[ $string67 \] $PAR +basis(1,7) $BR +\[ $string17 \] $PAR +basis(1,2) $BR +\[ $string12 \] $PAR +basis(2,3) $BR +\[ $string23\] $PAR +basis(3,2) $BR +\[ $string32 \] $PAR + +update tableau -- evaluation within BEGIN_TEXT/END_TEXT $PAR +basis(6,7)$BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(6,7)) )\} \] $PAR +basis(1,7) $BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(1,7)) )\} \] $PAR +basis(1,2) $BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(1,2)) )\} \] $PAR +basis(2,3) $BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(2,3)) )\}\] $PAR +basis(3,2) $BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(3,2)) )\} \] $PAR + + + + +END_TEXT +Context()->normalStrings; +TEXT($END_ONE_COLUMN); +ENDDOCUMENT(); diff --git a/t/matrix_tableau_tests/tableau-test6c.pg b/t/matrix_tableau_tests/tableau-test6c.pg new file mode 100644 index 0000000000..6db437b1ec --- /dev/null +++ b/t/matrix_tableau_tests/tableau-test6c.pg @@ -0,0 +1,237 @@ +DOCUMENT(); +loadMacros( +"PGstandard.pl", +"MathObjects.pl", +"parserPopUp.pl", +"unionLists.pl", +"MatrixReduce.pl", +#"AppletObjects.pl", +"PGessaymacros.pl", +"PGmatrixmacros.pl", +"PGML.pl", +"LinearProgramming.pl", +"parserLinearInequality.pl", +"quickMatrixEntry.pl", +#"scaffold.pl", +"tableau.pl", +#"gage_matrix_ops.pl", +"PGinfo.pl", +"source.pl", +"PGcourse.pl", +); + +TEXT(beginproblem()); +TEXT($BEGIN_ONE_COLUMN); +$showPartialCorrectAnswers = 1; + +INITIALIZE_QUICK_MATRIX_ENTRY(); + +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + + +#Hack to prevent domain conflict in answer. +Context()->variables->add(p1=>'Real',p2=>'Real'); +$objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); +# why can't the formula be defined within context "linearInequality"? + +Context("LinearInequality"); +Context()->variables->add(p1=>'Real',p2=>'Real'); +Context()->strings->add("Essay Answer" =>{}); +Context()->strings->add('Minimize'=>{},'Maximize'=>{}, "?"=>{}); +Context()->strings->add('Yes'=>{},'No'=>{}); +Context()->flags->set({tolType=>"absolute",tolerance=>.001}); +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +# $objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); + +our $context=Context(); + +sub basis_size { + my $basis = shift; + return scalar( @{$basis->data} ); +} +sub basis_matrix { + my $original_matrix = shift; + my $basis = shift; + my $basis_size=basis_size($basis); + my $basis_matrix = $original_matrix->submatrix(rows=>[1..$basis_size],columns=>($basis->data)); + return $basis_matrix; +} + +sub update_tableau { + my $original_matrix = shift; + my $basis = shift; #basis columns -- a Set + # $basis is a Set object + # last row is numbered $basis_size+1; + my $context= Context(); + #Context()->normalStrings; + my $basis_matrix = basis_matrix($original_matrix, $basis); + my $basis_size=basis_size($basis); + my $det = Real($basis_matrix->det)->value; # prevent tex output from det + my $obj_coeff = - ($original_matrix->row($basis_size+1)); + #obj function coeff corresponding to basis variables/columns + $cB = Matrix($obj_coeff->column_slice($basis->data)); + # get A | S |b portion of tableau and multiply by det(B)B^-1 + my $current_matrix= ($det)*($basis_matrix->inverse)*($original_matrix->row_slice(1..$basis_size)); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + #Context()->texStrings; # restore context + return $current_tableau; +} +# the full tableau has one more row -- the objective function +sub display_tableau_state { + my ($original_matrix,$tableau, $matrix, $basis,$pivot)= @_; + $basis = $basis->sort; + my $basis_size = @{$basis->data} ; + my $basis_matrix = matrix_from_submatrix($original_matrix, rows=>[1..$basis_size],columns=>$basis->data); + my $normalized_tableau = $matrix; + my $reduced_matrix = matrix_from_matrix_rows($matrix,1..$basis_size); + warn "basis matrix is $basis_matrix, $reduced_matrix"; + my $normalized_reduced_matrix = ($basis_matrix->det)*$reduced_matrix; + # calculate current tableau by multiplying by the inverse of the basis + my $obj_coeff = - $original_matrix->row($basis_size+1); + my $cB = Matrix($obj_coeff->column_slice($basis->data)); + my $current_matrix = ($basis_matrix->det)*($basis_matrix->inverse)*matrix_from_matrix_rows($original_matrix,1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + $current_tableau = update_tableau($original_matrix, $basis); + + #DEBUG_MESSAGE("get state variables for $matrix using basis $basis"); + @statevars1 = get_tableau_variable_values($matrix, $basis); + # get z value + $statevars1 = ~~@statevars1; + $state = Matrix([[@statevars1]]); + + + return " + + pivot: \($pivot\) basis: \( $basis\) $PAR basis matrix \(". + display_matrix_mm($basis_matrix)."\) $PAR + tableau and normalized current tableau $PAR + \(" . lp_display_mm($tableau) . "\)\(".lp_display_mm($normalized_tableau). " \) $PAR + original matrix and current matrix calculated by Binverse*original $PAR + \(". + display_matrix_mm(matrix_from_matrix_rows($original_matrix,1..$basis_size)). + "\) \(". + display_matrix_mm($current_matrix). + "\) $PAR + objective coefficients = \( $obj_coeff \) $PAR + basis objective_coefficients = \($cB\) $PAR + new tableau \[". + lp_display_mm($current_tableau). + "\] $PAR + state: \($state)\) $PAR + " +} + + +#### problem starts here: + +Context($context); +#constraint matrix +$A = Matrix([[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0], + [ -1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, ]]); +$b = Matrix([[-$bill_profit,-$steve_profit]])->transpose; +$c = Matrix([1,0,0,0,0]); + +my $tableau1 = Tableau->new(A=>$A, b=>$b, c=>$c); +$m = $tableau1->{M}; +$k = $tableau1->current_tableau; +$basis_cols24 = $tableau1->basis(2,4); + +$k24 = $tableau1->current_tableau; +$basis_cols42 = $tableau1->basis(4,2); + +$k42 = $tableau1->current_tableau; + +Context()->texStrings; +BEGIN_PGML + +A [`[$A]`], b [`[$b]`], C [`[$c]`] + +dimensions = [` [@ join(" ", $b->row(1)) @] `] +[` [@ join(" ", $b->column(1)) @] `] + +tableau [` [$tableau1->{M}] `] + +tableau m [`[$m]`] + +tableau k [`[$k]`] + +---------------------------------- +basis columns [$basis_cols24] + +tableau k24 [`[$k24]`] + +basis columns [$basis_cols42] + +tableau k42 [`[$k42]`] + +END_PGML +Context()->normalStrings; + +ENDDOCUMENT + + + +$ra_matrix = [[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1,0,0,0,-$bill_profit], + [-1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, 0,1,0,0,-$steve_profit], + [1,0,0,0,0,0,0,0,1,0]]; + +# [0,$money_total,$time_total,1,1,0,0,1,0,0]]; +$original_matrix = Matrix($ra_matrix); + +$ra_matrix1 = lp_clone($ra_matrix); + + +$string67 = lp_display_mm( update_tableau($original_matrix, Set(6,7)) ); +$string17 = lp_display_mm( update_tableau($original_matrix, Set(1,7)) ); +$string12 = lp_display_mm( update_tableau($original_matrix, Set(1,2)) ); +$string23 = lp_display_mm( update_tableau($original_matrix, Set(2,3)) ); +$string32 = lp_display_mm( update_tableau($original_matrix, Set(3,2)) ); + +Context()->texStrings; +BEGIN_TEXT +update tableau method tableau-test6b$PAR +basis(6,7)$BR +\[ $string67 \] $PAR +basis(1,7) $BR +\[ $string17 \] $PAR +basis(1,2) $BR +\[ $string12 \] $PAR +basis(2,3) $BR +\[ $string23\] $PAR +basis(3,2) $BR +\[ $string32 \] $PAR + +update tableau -- evaluation within BEGIN_TEXT/END_TEXT $PAR +basis(6,7)$BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(6,7)) )\} \] $PAR +basis(1,7) $BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(1,7)) )\} \] $PAR +basis(1,2) $BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(1,2)) )\} \] $PAR +basis(2,3) $BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(2,3)) )\}\] $PAR +basis(3,2) $BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(3,2)) )\} \] $PAR + + + + +END_TEXT +Context()->normalStrings; +TEXT($END_ONE_COLUMN); +ENDDOCUMENT(); diff --git a/t/matrix_tableau_tests/tableau-test6d.pg b/t/matrix_tableau_tests/tableau-test6d.pg new file mode 100644 index 0000000000..37f3b4bc43 --- /dev/null +++ b/t/matrix_tableau_tests/tableau-test6d.pg @@ -0,0 +1,206 @@ +DOCUMENT(); +loadMacros( +"PGstandard.pl", +"MathObjects.pl", +"parserPopUp.pl", +"unionLists.pl", +"MatrixReduce.pl", +#"AppletObjects.pl", +"PGessaymacros.pl", +"PGmatrixmacros.pl", +"LinearProgramming.pl", +"parserLinearInequality.pl", +"quickMatrixEntry.pl", +#"scaffold.pl", +"tableau.pl", +#"gage_matrix_ops.pl", +"PGinfo.pl", +"source.pl", +"PGcourse.pl", +); + +TEXT(beginproblem()); +TEXT($BEGIN_ONE_COLUMN); +$showPartialCorrectAnswers = 1; + +INITIALIZE_QUICK_MATRIX_ENTRY(); + +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + + +#Hack to prevent domain conflict in answer. +Context()->variables->add(p1=>'Real',p2=>'Real'); +$objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); +# why can't the formula be defined within context "linearInequality"? + +Context("LinearInequality"); +Context()->variables->add(p1=>'Real',p2=>'Real'); +Context()->strings->add("Essay Answer" =>{}); +Context()->strings->add('Minimize'=>{},'Maximize'=>{}, "?"=>{}); +Context()->strings->add('Yes'=>{},'No'=>{}); +Context()->flags->set({tolType=>"absolute",tolerance=>.001}); +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +# $objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); + +our $context=Context(); + +sub basis_size { + my $basis = shift; + return scalar( @{$basis->data} ); +} +sub basis_matrix { + my $original_matrix = shift; + my $basis = shift; + my $basis_size=basis_size($basis); + my $basis_matrix = $original_matrix->submatrix(rows=>[1..$basis_size],columns=>($basis->data)); + return $basis_matrix; +} + +sub update_tableau { + my $original_matrix = shift; + my $basis = shift; #basis columns -- a Set + # $basis is a Set object + # last row is numbered $basis_size+1; + my $context= Context(); + #Context()->normalStrings; + my $basis_matrix = basis_matrix($original_matrix, $basis); + my $basis_size=basis_size($basis); + my $det = Real($basis_matrix->det)->value; # prevent tex output from det + my $obj_coeff = - ($original_matrix->row($basis_size+1)); + #obj function coeff corresponding to basis variables/columns + $cB = Matrix($obj_coeff->column_slice($basis->data)); + # get A | S |b portion of tableau and multiply by det(B)B^-1 + my $current_matrix= ($det)*($basis_matrix->inverse)*($original_matrix->row_slice(1..$basis_size)); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + #Context()->texStrings; # restore context + return $current_tableau; +} +# the full tableau has one more row -- the objective function +sub display_tableau_state { + my ($original_matrix,$tableau, $matrix, $basis,$pivot)= @_; + $basis = $basis->sort; + my $basis_size = @{$basis->data} ; + my $basis_matrix = matrix_from_submatrix($original_matrix, rows=>[1..$basis_size],columns=>$basis->data); + my $normalized_tableau = $matrix; + my $reduced_matrix = matrix_from_matrix_rows($matrix,1..$basis_size); + warn "basis matrix is $basis_matrix, $reduced_matrix"; + my $normalized_reduced_matrix = ($basis_matrix->det)*$reduced_matrix; + # calculate current tableau by multiplying by the inverse of the basis + my $obj_coeff = - $original_matrix->row($basis_size+1); + my $cB = Matrix($obj_coeff->column_slice($basis->data)); + my $current_matrix = ($basis_matrix->det)*($basis_matrix->inverse)*matrix_from_matrix_rows($original_matrix,1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + $current_tableau = update_tableau($original_matrix, $basis); + + #DEBUG_MESSAGE("get state variables for $matrix using basis $basis"); + @statevars1 = get_tableau_variable_values($matrix, $basis); + # get z value + $statevars1 = ~~@statevars1; + $state = Matrix([[@statevars1]]); + + + return " + + pivot: \($pivot\) basis: \( $basis\) $PAR basis matrix \(". + display_matrix_mm($basis_matrix)."\) $PAR + tableau and normalized current tableau $PAR + \(" . lp_display_mm($tableau) . "\)\(".lp_display_mm($normalized_tableau). " \) $PAR + original matrix and current matrix calculated by Binverse*original $PAR + \(". + display_matrix_mm(matrix_from_matrix_rows($original_matrix,1..$basis_size)). + "\) \(". + display_matrix_mm($current_matrix). + "\) $PAR + objective coefficients = \( $obj_coeff \) $PAR + basis objective_coefficients = \($cB\) $PAR + new tableau \[". + lp_display_mm($current_tableau). + "\] $PAR + state: \($state)\) $PAR + " +} + + +#### problem starts here: + +Context($context); + +$ra_matrix = [[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1,0,0,0,-$bill_profit], + [-1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, 0,1,0,0,-$steve_profit], + [1,0,0,0,0,0,0,0,1,0]]; + +# [0,$money_total,$time_total,1,1,0,0,1,0,0]]; +$original_matrix = Matrix($ra_matrix); + +$ra_matrix1 = lp_clone($ra_matrix); + +$A = Matrix([[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0], + [ -1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, ]]); +$b = Matrix([[-$bill_profit,-$steve_profit]])->transpose; +$c = Matrix([1,0,0,0,0]); + +$tableau1 = Tableau->new(A=>$A, b=>$b, c=>$c); + +# [0,$money_total,$time_total,1,1,0,0,1,0,0]]; + + + +$basis_cols67 = $tableau1->basis(6,7); +$string67 = lp_display_mm( $tableau1->current_tableau ); +$tableau1->basis(1,7); +$string17 = lp_display_mm( $tableau1->current_tableau ); +$tableau1->basis(1,2); +$string12 = lp_display_mm( $tableau1->current_tableau ); +$tableau1->basis(2,3); +$string23 = lp_display_mm( $tableau1->current_tableau ); +$tableau1->basis(3,2); +$string32 = lp_display_mm( $tableau1->current_tableau ); + +Context()->texStrings; +BEGIN_TEXT +update tableau method tableau-test6b$PAR +basis(6,7)$BR +\[ $string67 \] $PAR +basis(1,7) $BR +\[ $string17 \] $PAR +basis(1,2) $BR +\[ $string12 \] $PAR +basis(2,3) $BR +\[ $string23\] $PAR +basis(3,2) $BR +\[ $string32 \] $PAR + +update tableau -- evaluation within BEGIN_TEXT/END_TEXT $PAR +basis(6,7)$BR +\[ \{$tableau1->basis(6,7),lp_display_mm( $tableau1->current_tableau)\} \] $PAR +basis(1,7) $BR +\[ \{$tableau1->basis(1,7),lp_display_mm( $tableau1->current_tableau)\} \] $PAR +basis(1,2) $BR +\[ \{$tableau1->basis(1,2),lp_display_mm( $tableau1->current_tableau)\} \] $PAR +basis(2,3) $BR +\[ \{ $tableau1->basis(2,3),lp_display_mm( $tableau1->current_tableau)\} \] $PAR +basis(3,2) $BR +\[ \{lp_display_mm( update_tableau($original_matrix, Set(3,2)) )\} \] $PAR + + + + +END_TEXT +Context()->normalStrings; +TEXT($END_ONE_COLUMN); +ENDDOCUMENT(); diff --git a/t/matrix_tableau_tests/tableau-test6e.pg b/t/matrix_tableau_tests/tableau-test6e.pg new file mode 100644 index 0000000000..1b86699989 --- /dev/null +++ b/t/matrix_tableau_tests/tableau-test6e.pg @@ -0,0 +1,109 @@ +DOCUMENT(); +loadMacros( +"PGstandard.pl", +"MathObjects.pl", +"PGmatrixmacros.pl", +"niceTables.pl", +"tableau.pl", +"PGinfo.pl", +"source.pl", +"PGcourse.pl", +); + +TEXT(beginproblem()); +TEXT($BEGIN_ONE_COLUMN); +$showPartialCorrectAnswers = 1; + + +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + + + +#### problem starts here: + +# need error checking to make sure that tableau->new checks +# that inputs are matrices + +$A = Matrix([[-$bill_money_commitment,-$bill_time_commitment, -1, 0], + [ -$steve_money_commitment,-$steve_time_commitment, 0, -1 ]]); +$b = Vector([-$bill_profit,-$steve_profit]); # need vertical vector +$c = ColumnVector([$money_total,$time_total,0,0]); + +BEGIN_TEXT +A => $A $PAR +b => $b $PAR +c => $c $PAR +END_TEXT +# $c = Matrix([[$money_total,$time_total,0,0],[-1,0,0,0,0]]); +# not ready for multiple objective rows + +$tableau1 = Tableau->new(A=>$A, b=>$b, c=>$c); +# slack variables are automatically added + + + +$toprow = [qw(y1 y2 y3 y4 s1 s2 P b) ]; +$align = "cccc|cc|c|c"; + +BEGIN_TEXT +tableau is $BR \(\{$tableau1\}\)$BR +tableau->current_tableau is $BR +\(\{$tableau1->current_tableau\}\)$BR +tableau {current_constraint_matrix} $BR \(\{$tableau1->{current_constraint_matrix}\}\)$BR +tableau {M} $BR \(\{$tableau1->{M}\}\) $BR +current objective row coefficients of tableau $BR \(\{$tableau1->{current_coeff}\}\)$PAR +lop_display() pretty prints the value of tableau->current_tableau with optional decorations. +The usual decorations are an alignment scheme which allows one to add vertical lines +and labels for the rows of the matrix (on the top). The last row is automatically separated +off as the line storing the (negative of) the objective coefficients. + +lop_display(tableau) with no options set $BR +\{lop_display($tableau1)\}$BR +lop_display(tableau, alignment=>"$align", toplevel=>[qw(y1 y2 y3 y4 s1 s2 P b) ])$BR +\{lop_display($tableau1,alignment=>$align,toplevel=>$toprow)\}$PAR + +ok? +$HR +$PAR + +$PAR +END_TEXT + +BEGIN_TEXT +$PAR +OTHER RANDOM TESTS +$PAR + +obj row \{$tableau1->{current_coeff}\} +$PAR +tableau looks like \{pretty_print( $tableau1 )\} + + + + + +$PAR +with datatable +$PAR +\{DataTable([@matrix,[[3,midrule=>1],2,0,0,0,0,1,0]], caption=>"Values",align=>"ccccc|cc|c|c") \} + +$PAR +\{LayoutTable([@matrix,[[3,midrule=>1],2,0,0,0,0,1,0]], caption=>"Values",align=>"ccccc|cc|c|c") \} + +$PAR +\{DataTable( [[1,2,3],[4,5,6]], caption=>"Values") \} + +END_TEXT +Context()->normalStrings; +TEXT($END_ONE_COLUMN); +ENDDOCUMENT(); diff --git a/t/matrix_tableau_tests/tableau_answer_test1.pg b/t/matrix_tableau_tests/tableau_answer_test1.pg new file mode 100644 index 0000000000..753b46cc12 --- /dev/null +++ b/t/matrix_tableau_tests/tableau_answer_test1.pg @@ -0,0 +1,105 @@ +##DESCRIPTION + + + +##ENDDESCRIPTION + + +DOCUMENT(); # This should be the first executable line in the problem. + +loadMacros( + "PGstandard.pl", # Standard macros for PG language + "MathObjects.pl", + "PGML.pl", + "tableau.pl", + #"source.pl", # used to display problem source button + "PGcourse.pl", # Customization file for the course +); + +TEXT(beginproblem()); +$showPartialCorrectAnswers = 1; + +############################################################## +# +# Setup +# +# +Context("Matrix"); + +$m = Matrix([1,2,3],[4,5,6],[7,8,9]); +$n = Matrix([4,5,6],[2,4,6],[7,8,9]); + +@rows1 = matrix_extract_rows($m); +@rows2 = matrix_extract_rows($n); + +@rows1a= map {Vector($_)} @rows1; +TEXT( join(" ", @rows1a)); +@rows2a= map {Vector($_)} @rows2; +############################################################## +# +# Text +# +# +sub tableauEquivalence_orig { + return sub { + my ($correct,$student,$ansHash,$nth,$value)=@_; + my $score = Vector($correct)->cmp(parallel=>1)->evaluate(Vector($student))->{score}; + return $score; + } + } + +# $rh_ans=List(@rows1a)->cmp(checker=>sub { +# my ($correct,$student,$ansHash,$nth,$value)=@_; +# my $score = $correct->cmp(parallel=>1)->evaluate($student)->{score}; +# return $score; +# })->evaluate(List(@rows2a)); + +$rh_ans=List(@rows1a)->cmp(checker=>tableauEquivalence_orig())->evaluate(List(@rows2a)); +$rh_ans2 = List(@rows1)->cmp(checker=>tableauEquivalence_orig())->evaluate(List(@rows2)); +Context()->texStrings; + +$rh_ans3 = $m->cmp(checker=>tableauEquivalence())->evaluate($n); +BEGIN_TEXT +rows : \[\{join(" ", @rows1,"|||",@rows2)\}\]; $PAR + +rows are equal \{$rh_ans->{score}\} $PAR +rows are equal (even with out being classed to Vector) \{$rh_ans2->{score}\} $PAR +Comparing as matrices \{$rh_ans3->{score}\}$PAR +\{pretty_print($rh_ans)\} $PAR +\{pretty_print($rh_ans2)\}$PAR + + + +$PAR + +END_TEXT + +Context()->normalStrings; + +# This works in online PG lab + +# $a= Vector(1,2,3); +# $b= Vector(2,4,6); +# $c= Vector(4,2,3); +# $d=Vector(8,4,6); +# +# $list2=List($d,$a); +# TEXT( $list1 = List($a,$c)->cmp( +# checker=> sub { +# my ($correct,$student,$ansHash,$nth,$value)=@_; +# my $score = $correct->cmp(parallel=>1)->evaluate($student)->{score}; +# return $score; +# } +# )->evaluate($list2)->pretty_print ); +# +# TEXT($a->cmp(parallel=>1)->evaluate($b)->pretty_print); + +############################################################## +# +# Answers +# +# + + + +ENDDOCUMENT(); # This should be the last executable line in the problem. \ No newline at end of file diff --git a/t/matrix_tableau_tests/tableau_module_test1.pg b/t/matrix_tableau_tests/tableau_module_test1.pg new file mode 100644 index 0000000000..30cf72eb2f --- /dev/null +++ b/t/matrix_tableau_tests/tableau_module_test1.pg @@ -0,0 +1,79 @@ +############################################## +DOCUMENT(); + +loadMacros( + "PGstandard.pl", # Standard macros for PG language + "MathObjects.pl", + "parserLinearInequality.pl", + "PGML.pl", + "tableau.pl", + "PGmatrixmacros.pl", + "LinearProgramming.pl", + #"source.pl", # allows code to be displayed on certain sites. + "PGcourse.pl", +); + +############################################## + +Context("Matrix"); # need Matrix context to allow string input into Matrix. + +$m = Matrix("[[3,6,7],[2,1,8],[4,6,21],[-6,7,9]]"); +$constraint_matrix = Matrix(" +[[ 0, 0, -1, -1], + [-1, -1, 0, 0 ], + [1, 0 , 1 , 0], + [0, 1, 0, 1]] +"); + +#TEXT ("created ". ref($m)); +#what are the best ways to display a matrix? + +$m1 = $m->row_slice([4,1]); +$m2 = $m->column_slice([3,2,1]); + +$list = $m->extract_rows_to_list(2,3); + +$b = Matrix([1, 2, 3, 4])->transpose; # make it an n by 1 matrix +$c = Matrix([5, 6, 7]); +$t = Tableau->new(A=>$m,b=>$b, c=>$c); + +$c_4 = $t->current_tableau; +############################################################## +# +# Text +# +# + +Context()->texStrings; + +# tableau is \[$c_4\] $PAR +# +# inline tableau is \[\{ $t->current_tableau \}\] +# +# tableau->{M} is \[\{ $t->{M} \}\] +# +# $PAR +# Enter tableau->{M} \{$t->{M}->ans_array\} + +BEGIN_TEXT +m is \[$m\] + +$PAR +enter m = \{$m->ans_array\} + +\{#$m->cmp->evaluate($m)->pretty_print\} + +END_TEXT +Context()->normalStrings; + + +############################################################## +# +# Answers +# +# +#ANS($t->{M}->cmp(checker=>tableauEquivalence()) ); +ANS($m->cmp(checker=>tableauEquivalence())); + + +ENDDOCUMENT(); # This should be the last executable line in the problem. \ No newline at end of file diff --git a/t/matrix_tableau_tests/tableau_row_operations_test.pg b/t/matrix_tableau_tests/tableau_row_operations_test.pg new file mode 100644 index 0000000000..c449ea1838 --- /dev/null +++ b/t/matrix_tableau_tests/tableau_row_operations_test.pg @@ -0,0 +1,206 @@ +##DESCRIPTION + + + +##ENDDESCRIPTION + + +DOCUMENT(); # This should be the first executable line in the problem. + +loadMacros( + "PGstandard.pl", # Standard macros for PG language + "MathObjects.pl", + "MatrixReduce.pl", + #"AppletObjects.pl", + "PGessaymacros.pl", + "PGmatrixmacros.pl", + "PGML.pl", + "LinearProgramming.pl", + "parserLinearInequality.pl", + "quickMatrixEntry.pl", + #"scaffold.pl", + "tableau.pl", + #"gage_matrix_ops.pl", + "PGinfo.pl", + "source.pl", + "PGcourse.pl", +); + +TEXT(beginproblem()); +$showPartialCorrectAnswers = 1; + +############################################################## +# +# Setup +# +# +Context("Matrix"); + +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +INITIALIZE_QUICK_MATRIX_ENTRY(); + +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + + +#constraint matrix +$A = Matrix([[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0], + [ -1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, ]]); +$b = Matrix([[-$bill_profit,-$steve_profit]])->transpose; +$c = Matrix([1,0,0,0,0]); + +$tableau1 = Tableau->new(A=>$A, b=>$b, c=>$c); +$m = $tableau1->{M}; +$m2 = $tableau1->{current_constraint_matrix}; +$tableau_orig = $tableau1->current_tableau; +$basis_orig = $tableau1->basis; +#HR +$tableau1->basis(1,2); +$m3 = $tableau1->{M}; +$m4 = $tableau1->{current_constraint_matrix}; +$tableau_orig3 = $tableau1->current_tableau; +$basis_orig3 = $tableau1->basis; +#HR +$tableau1->basis(6,7); +$t1 = $tableau1->current_tableau; +($index1, $value1, $feasible1) = $tableau1->find_short_cut_row(); +($index2, $value2, $infeasible2) = $tableau1->find_short_cut_column(1); +# search row 1 +($index3, $value3) = $tableau1->find_leaving_column(1); +$basis17 = $tableau1->find_next_basis_from_pivot(1,1); +$basis17=Set($basis17)->sort ; + +$tableau1->basis($basis17->value); # basis should be able to accept Set and List objects +$basis17check = $tableau1->basis; +$tableau17= $tableau1->current_tableau; +$matrix17 = $tableau1->{current_constraint_matrix}; # M = A | S |b + +#HR last part of phase 1 +($index11, $value11, $feasible11) = $tableau1->find_short_cut_row(); + +#HR first pivot of phase2 +($pivot1col, $value1c, $optimum1a) = $tableau1->find_pivot_column('min'); +($pivot1row, $value1r, $unbounded1a) = $tableau1->find_pivot_row(2); +#($row2,$col2,$optimum2,$unbounded2) = $tableau1->find_next_pivot('min'); +#$basis2 = $tableau1->find_next_basis_from_pivot($row2, $col2); +#TEXT("basis2 has type ".ref($basis2)); +#$tableau1->basis($basis2->value); +#$tableaudisplay1 = $tableau1->current_tableau; + +#HR second pivot of phase2 +#($pivot3col, $value3c, $optimum3a) = $tableau1->find_pivot_column('min'); +#($pivot3row, $value3r, $unbounded3a) = $tableau1->find_pivot_row(3); +#($row3,$col3,$optimum3,$unbounded3) = $tableau1->find_next_pivot('min'); +#$basis3 = $tableau1->find_next_basis_from_pivot($row3, $col3); +#$tableau1->basis($basis3->value); +#$tableaudisplay2 = $tableau1->current_tableau; + +#HR third pivot of phase2 +#($pivot4col, $value4c, $optimum4a) = $tableau1->find_pivot_column('min'); +#($pivot4row, $value4r, $unbounded4a) = $tableau1->find_pivot_row(3); +#($row4,$col4,$optimum4,$unbounded4) = $tableau1->find_next_pivot('min'); +#optimum is reached so can't go forward from here +#$basis4 = $tableau1->find_next_basis_from_pivot($row4, $col4); +#$tableau1->basis($basis4->value); +#$tableaudisplay3 = $tableau1->current_tableau; + +### experimental +# basis calculation +TEXT(" start with basis (6,7) and pivot on (1,1)"); + $basis = Set(6,7); + $basis = $basis - Set(6); + $basis = List( 1, $basis->value, ); +TEXT( "new basis is $basis"); + +############################################################## +# +# Text +# +# + +Context()->texStrings; +BEGIN_TEXT +matrix: \[ $m \quad $m2 \] $PAR +tableau: \[ $tableau_orig \] $PAR +original basis: $basis_orig $PAR +$HR +matrix: \[ $m3 \quad $m4 \] $PAR +tableau: \[ $tableau_orig3 \] $PAR +new basis: $basis_orig3 $PAR +$HR +\[ $t1\] +find shortcut row: (index, value, feasible) = ($index1, $value1, $feasible1) $PAR +Since it is not feasible yet we'll search row 1 for a short cut pivot $PAR +find shortcut column: (index, value, infeasible) = ($index2, $value2, $infeasible2) $PAR +find leaving_column for row 1: (col_index,value) = ($index3, $value3) $PAR +find new basis from pivot at (1,1): \($basis17\) $PAR +check basis: ($basis17check), \($basis17check\) $PAR +tableau17, matrix17: \[ $tableau17 \qquad $matrix17 \] $PAR +$HR +last part of phase 1 $PAR +find shortcut row: (row $index11,value $value11,feasible $feasible11) $PAR +since the tableau is now feasible we move on to phase2 +$HR + +first pivot of phase 2: $PAR +find pivot col ($pivot1col, $value1c, $optimum1) $PAR +find pivot row ($pivot1row, $value1r, $unbounded1) $PAR +next pivot:(row col optimum unbounded)= +( \{join(",",$tableau1->find_next_pivot('min'))\} ) + $PAR +next basis: (\{$tableau1->basis( +$tableau1->find_next_basis_from_pivot( +($tableau1->find_next_pivot('min'))[0], +($tableau1->find_next_pivot('min'))[1] +)->value +)\}) $PAR +new tableau: $PAR +\[ \{ $tableau1->current_tableau \} \]$PAR + +$HR + +second pivot of phase2: $PAR +next pivot:(row col optimum unbounded)= +( \{join(",",$tableau1->find_next_pivot('min'))\} ) + $PAR +next basis: (\{$tableau1->basis( +$tableau1->find_next_basis_from_pivot( +($tableau1->find_next_pivot('min'))[0], +($tableau1->find_next_pivot('min'))[1] +)->value +)\}) $PAR +new tableau: $PAR +\[ \{ $tableau1->current_tableau \} \]$PAR + +$HR + +third pivot of phase2: $PAR +next pivot:(row col optimum unbounded)= +( \{join(",",$tableau1->find_next_pivot('min'))\} ) +optimum so we stop here + +END_TEXT +Context()->normalStrings; + + +############################################################## +# +# Answers +# +# + + + +ENDDOCUMENT(); # This should be the last executable line in the problem. \ No newline at end of file diff --git a/t/matrix_tableau_tests/tableau_row_operations_testB.pg b/t/matrix_tableau_tests/tableau_row_operations_testB.pg new file mode 100644 index 0000000000..6d681714e9 --- /dev/null +++ b/t/matrix_tableau_tests/tableau_row_operations_testB.pg @@ -0,0 +1,325 @@ +##DESCRIPTION + + + +##ENDDESCRIPTION + + +DOCUMENT(); # This should be the first executable line in the problem. + +loadMacros( + "PGstandard.pl", # Standard macros for PG language + "MathObjects.pl", + "MatrixReduce.pl", + #"AppletObjects.pl", + "PGessaymacros.pl", + "PGmatrixmacros.pl", + "PGML.pl", + "LinearProgramming.pl", + "parserLinearInequality.pl", + "quickMatrixEntry.pl", + #"scaffold.pl", + "tableau.pl", + #"gage_matrix_ops.pl", + "PGinfo.pl", + "source.pl", + "PGcourse.pl", +); + +TEXT(beginproblem()); +$showPartialCorrectAnswers = 1; + +############################################################## +# +# Setup +# +# +Context("Matrix"); + +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +INITIALIZE_QUICK_MATRIX_ENTRY(); + +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + + +#constraint matrix +$A = Matrix([[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0], + [ -1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, ]]); +$b = Matrix([[-$bill_profit,-$steve_profit]])->transpose; +$c = Matrix([1,0,0,0,0]); + +$tableau1 = Tableau->new(A=>$A, b=>$b, c=>$c); +$m = $tableau1->{M}; +$m2 = $tableau1->{current_constraint_matrix}; +$tableau_orig = $tableau1->current_tableau; +$basis_orig = $tableau1->basis; +#HR +$tableau1->basis(1,2); +$m3 = $tableau1->{M}; +$m4 = $tableau1->{current_constraint_matrix}; +$tableau_orig3 = $tableau1->current_tableau; +$basis_orig3 = $tableau1->basis; +#HR +$tableau1->basis(6,7); +$t1 = $tableau1->current_tableau; +($index1, $value1, $feasible1) = $tableau1->find_short_cut_row(); +($index2, $value2, $infeasible2) = $tableau1->find_short_cut_column(1); +# search row 1 +($index3, $value3) = $tableau1->find_leaving_column(1); +$basis17 = $tableau1->find_next_basis_from_pivot(1,1); +$basis17=Set($basis17)->sort ; + +$tableau1->basis($basis17->value); # basis should be able to accept Set and List objects +$basis17check = $tableau1->basis; +$tableau17= $tableau1->current_tableau; +$matrix17 = $tableau1->{current_constraint_matrix}; # M = A | S |b +#HR last part of phase 1 +($index11, $value11, $feasible11) = $tableau1->find_short_cut_row(); + +#HR first pivot of phase2 +($pivot1col, $value1c, $optimum1a) = $tableau1->find_pivot_column('min'); +($pivot1row, $value1r, $unbounded1a) = $tableau1->find_pivot_row(2); +#($row2,$col2,$optimum2,$unbounded2) = $tableau1->find_next_pivot('min'); +#$basis2 = $tableau1->find_next_basis_from_pivot($row2, $col2); +#TEXT("basis2 has type ".ref($basis2)); +#$tableau1->basis($basis2->value); +#$tableaudisplay1 = $tableau1->current_tableau; + +#HR second pivot of phase2 +#($pivot3col, $value3c, $optimum3a) = $tableau1->find_pivot_column('min'); +#($pivot3row, $value3r, $unbounded3a) = $tableau1->find_pivot_row(3); +#($row3,$col3,$optimum3,$unbounded3) = $tableau1->find_next_pivot('min'); +#$basis3 = $tableau1->find_next_basis_from_pivot($row3, $col3); +#$tableau1->basis($basis3->value); +#$tableaudisplay2 = $tableau1->current_tableau; + +#HR third pivot of phase2 +#($pivot4col, $value4c, $optimum4a) = $tableau1->find_pivot_column('min'); +#($pivot4row, $value4r, $unbounded4a) = $tableau1->find_pivot_row(3); +#($row4,$col4,$optimum4,$unbounded4) = $tableau1->find_next_pivot('min'); +#optimum is reached so can't go forward from here +#$basis4 = $tableau1->find_next_basis_from_pivot($row4, $col4); +#$tableau1->basis($basis4->value); +#$tableaudisplay3 = $tableau1->current_tableau; + +### experimental +# basis calculation +TEXT(" start with basis (6,7) and pivot on (1,1)"); + $basis = Set(6,7); + $basis = $basis - Set(6); + $basis = List( 1, $basis->value, ); +TEXT( "new basis is $basis"); + +############################################################## +# +# Text +# +# +$tableau1->basis(6,7); +@output_find_next_basis=@{$tableau1->find_next_basis('min')}; + +Context()->texStrings; +BEGIN_TEXT +matrix: \[ $m \quad $m2 \] $PAR +tableau: \[ $tableau_orig \] $PAR +original basis: $basis_orig $PAR +$HR +new basis (3,2) +\[\{$tableau1->current_tableau(3,2)\}\] +basis \{$tableau1->basis()\} ,\{List($tableau1->{basis})\} $PAR +basis matrix \[\{$tableau1->{current_constraint_matrix}-> + submatrix(rows=>[1..($tableau1->{m})],columns=>$tableau1->{basis})\}\]; +$PAR +new basis (2,3) +\[\{$tableau1->current_tableau(2,3)\}\] +basis \{$tableau1->basis()\}, \{List($tableau1->{basis})\}$PAR +basis matrix \[\{$tableau1->{current_constraint_matrix}->submatrix(rows=>[1..($tableau1->{m})],columns=>$tableau1->{basis})\}\]; +$PAR +Should figure out why these don't reverse although internally they do. +$HR +Figure out outputs: $PAR +find_next_basis: @output_find_next_basis $PAR +test1 \{[$tableau1->find_next_basis('max')]->[0]\},$PAR +test2 \{[$tableau1->find_next_basis('max')]->[1] \},$PAR + + +$HR +next pivot 1 : \{ join(",",$tableau1->find_next_pivot('max'))\}$PAR +next basis \{ join(",",$tableau1->find_next_basis('max')->[0])\}$PAR +tableau: $PAR +\[\{$tableau1->current_tableau( + [$tableau1->find_next_basis('max')]->[0], + [$tableau1->find_next_basis('max')]->[1] +)\}\] +$HR +next pivot 2 : \{ join(",",$tableau1->find_next_pivot('max'))\}$PAR +next basis : \{ join(",",$tableau1->find_next_basis('max'))\}$PAR +tableau: $PAR +\[\{$tableau1->current_tableau( + [$tableau1->find_next_basis('max')]->[0], + [$tableau1->find_next_basis('max')]->[1] +)\}\] +$HR +next pivot : \(\{List($tableau1->find_next_pivot('max')) \}\)$PAR +next basis = \(\{List($tableau1->find_next_basis('max')) \}\)$PAR +tableau: $PAR +\[\{$tableau1->current_tableau( + [$tableau1->find_next_basis('max')]->[0], + [$tableau1->find_next_basis('max')]->[1] +)\}\] +$HR + +Unit tests $PAR + +next pivot column \{join(",",$tableau1->find_pivot_column('max'))\}$PAR +next pivot row \{ join(",",$tableau1->find_pivot_row(6))\}$PAR +next pivot \{ join(",",$tableau1->find_next_pivot('max'))\}$PAR +next basis \{ join(",",$tableau1->find_next_basis('max'))\}$PAR +test pivot \{ ($tableau1->find_next_pivot('max'))[2]\}$PAR + +LOP is unbounded $PAR +$HR +matrix: \[ $m3 \quad $m4 \] $PAR +tableau: \[ $tableau_orig3 \] $PAR +new basis: $basis_orig3 $PAR +$HR +\[ $t1\] +find shortcut row: (index, value, feasible) = ($index1, $value1, $feasible1) $PAR +Since it is not feasible yet we'll search row 1 for a short cut pivot $PAR +find shortcut column: (index, value, infeasible) = ($index2, $value2, $infeasible2) $PAR +find leaving_column for row 1: (col_index,value) = ($index3, $value3) $PAR +find new basis from pivot at (1,1): \($basis17\) $PAR +check basis: ($basis17check), \($basis17check\) $PAR +tableau17, matrix17: \[ $tableau17 \qquad $matrix17 \] $PAR +$HR +last part of phase 1 $PAR +find shortcut row: (row $index11,value $value11,feasible $feasible11) $PAR +since the tableau is now feasible we move on to phase2 +$HR + +first pivot of phase 2: $PAR +find pivot col ($pivot1col, $value1c, $optimum1) $PAR +find pivot row ($pivot1row, $value1r, $unbounded1) $PAR +next pivot:(row col optimum unbounded)= +( \{join(",",$tableau1->find_next_pivot('min'))\} ) + $PAR +next basis: (\{$tableau1->basis( +$tableau1->find_next_basis_from_pivot( +[$tableau1->find_next_pivot('min')]->[0], +[$tableau1->find_next_pivot('min')]->[1] +)->value +)\}) $PAR +new tableau: $PAR +\[ \{ $tableau1->current_tableau \} \]$PAR + +$HR + +second pivot of phase2: $PAR +next pivot:(row col optimum unbounded)= +( \{join(",",$tableau1->find_next_pivot('min'))\} ) + $PAR +next basis: (\{$tableau1->basis( +$tableau1->find_next_basis_from_pivot( +[$tableau1->find_next_pivot('min')]->[0], +[$tableau1->find_next_pivot('min')]->[1] +)->value +)\}) $PAR +new tableau: $PAR +\[ \{ $tableau1->current_tableau \} \]$PAR + +$HR + +third pivot of phase2: $PAR +next pivot:(row col optimum unbounded)= +( \{join(",",$tableau1->find_next_pivot('min'))\} ) +optimum so we stop here and start going toward the maximum: + + + +$HR + +$HR + + + +next pivot1 : \(\{List($tableau1->find_next_pivot('max')) \}\)$PAR +next basis = \(\{List($tableau1->find_next_basis('max')) \}\)$PAR +tableau: $PAR +\[\{$tableau1->current_tableau( + [$tableau1->find_next_basis('max')]->[0], + [$tableau1->find_next_basis('max')]->[1] +)\}\] +$HR +next pivot2 : \(\{List($tableau1->find_next_pivot('max')) \}\)$PAR +next basis = \(\{List($tableau1->find_next_basis('max')) \}\)$PAR +tableau: $PAR +\[\{$tableau1->current_tableau( + [$tableau1->find_next_basis('max')]->[0], + [$tableau1->find_next_basis('max')]->[1] +)\}\] +$HR +next pivot3 : \(\{List($tableau1->find_next_pivot('max')) \}\)$PAR +next basis = \(\{List($tableau1->find_next_basis('max')) \}\)$PAR +tableau: $PAR +\[\{$tableau1->current_tableau( + [$tableau1->find_next_basis('max')]->[0], + [$tableau1->find_next_basis('max')]->[1] +)\}\] +$HR +next pivot4 ; (\{ join(",", $tableau1->find_next_pivot('max') ) \} )$PAR +next pivot4 : \(\{List($tableau1->find_next_pivot('max')) \}\)$PAR +[~~$tableau1->find_next_pivot('max')]->[3] $PAR +element 4 unbounded : \(\{[$tableau1->find_next_pivot('max')]->[3] \}\)$PAR +element 3 optimal : \(\{[$tableau1->find_next_pivot('max')]->[2] \}\)$PAR +This means we should stop $PAR +next basis = \(\{List($tableau1->find_next_basis('max')) \}\)$PAR +tableau: $PAR +\[\{$tableau1->current_tableau( + [$tableau1->find_next_basis('max')]->[0], + [$tableau1->find_next_basis('max')]->[1] +)\}\] +$HR +next pivot : \(\{List($tableau1->find_next_pivot('max')) \}\)$PAR +next basis = \(\{List($tableau1->find_next_basis('max')) \}\)$PAR +tableau: $PAR +\[\{$tableau1->current_tableau( + [$tableau1->find_next_basis('max')]->[0], + [$tableau1->find_next_basis('max')]->[1] +)\}\] +$HR +next pivot : \(\{List($tableau1->find_next_pivot('max')) \}\)$PAR +next basis = \(\{List($tableau1->find_next_basis('max')) \}\)$PAR +tableau: $PAR +\[\{lp_display_mm( + $tableau1->current_tableau( + [$tableau1->find_next_basis('max')]->[0], + [$tableau1->find_next_basis('max')]->[1] + ) +)\}\] + +END_TEXT +Context()->normalStrings; + + +############################################################## +# +# Answers +# +# + + + +ENDDOCUMENT(); # This should be the last executable line in the problem. \ No newline at end of file diff --git a/t/matrix_tableau_tests/tableau_row_operations_testC.pg b/t/matrix_tableau_tests/tableau_row_operations_testC.pg new file mode 100644 index 0000000000..c207571f00 --- /dev/null +++ b/t/matrix_tableau_tests/tableau_row_operations_testC.pg @@ -0,0 +1,189 @@ +##DESCRIPTION + + + +##ENDDESCRIPTION + + +DOCUMENT(); # This should be the first executable line in the problem. + +loadMacros( + "PGstandard.pl", # Standard macros for PG language + "MathObjects.pl", + "MatrixReduce.pl", + #"AppletObjects.pl", + "PGessaymacros.pl", + "PGmatrixmacros.pl", + "PGML.pl", + "LinearProgramming.pl", + "parserLinearInequality.pl", + "quickMatrixEntry.pl", + #"scaffold.pl", + "tableau.pl", + #"gage_matrix_ops.pl", + "PGinfo.pl", + "source.pl", + "PGcourse.pl", +); + +TEXT(beginproblem()); +$showPartialCorrectAnswers = 1; + +############################################################## +# +# Setup +# +# +Context("Matrix"); + +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +INITIALIZE_QUICK_MATRIX_ENTRY(); + +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + + +#constraint matrix +$A = Matrix([[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0], + [ -1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, ]]); +$b = Matrix([[-$bill_profit,-$steve_profit]])->transpose; +$c = Matrix([1,0,0,0,0]); + +$tableau1 = Tableau->new(A=>$A, b=>$b, c=>$c); +$m = $tableau1->{M}; +$m2 = $tableau1->{current_constraint_matrix}; +$tableau_orig = $tableau1->current_tableau; +$basis_orig = $tableau1->basis; +#HR +$tableau1->basis(1,2); +$m3 = $tableau1->{M}; +$m4 = $tableau1->{current_constraint_matrix}; +$tableau_orig3 = $tableau1->current_tableau; +$basis_orig3 = $tableau1->basis; +#HR +$tableau1->basis(6,7); +$t1 = $tableau1->current_tableau; +($index1, $value1, $feasible1) = $tableau1->find_short_cut_row(); +($index2, $value2, $infeasible2) = $tableau1->find_short_cut_column(1); +# search row 1 +($index3, $value3) = $tableau1->find_leaving_column(1); +$basis17 = $tableau1->find_next_basis_from_pivot(1,1); +$basis17=Set($basis17)->sort ; + +$tableau1->basis($basis17->value); # basis should be able to accept Set and List objects +$basis17check = $tableau1->basis; +$tableau17= $tableau1->current_tableau; +$matrix17 = $tableau1->{current_constraint_matrix}; # M = A | S |b +#HR last part of phase 1 +($index11, $value11, $feasible11) = $tableau1->find_short_cut_row(); + +#HR first pivot of phase2 +($pivot1col, $value1c, $optimum1a) = $tableau1->find_pivot_column('min'); +($pivot1row, $value1r, $unbounded1a) = $tableau1->find_pivot_row(2); +#($row2,$col2,$optimum2,$unbounded2) = $tableau1->find_next_pivot('min'); +#$basis2 = $tableau1->find_next_basis_from_pivot($row2, $col2); +#TEXT("basis2 has type ".ref($basis2)); +#$tableau1->basis($basis2->value); +#$tableaudisplay1 = $tableau1->current_tableau; + +#HR second pivot of phase2 +#($pivot3col, $value3c, $optimum3a) = $tableau1->find_pivot_column('min'); +#($pivot3row, $value3r, $unbounded3a) = $tableau1->find_pivot_row(3); +#($row3,$col3,$optimum3,$unbounded3) = $tableau1->find_next_pivot('min'); +#$basis3 = $tableau1->find_next_basis_from_pivot($row3, $col3); +#$tableau1->basis($basis3->value); +#$tableaudisplay2 = $tableau1->current_tableau; + +#HR third pivot of phase2 +#($pivot4col, $value4c, $optimum4a) = $tableau1->find_pivot_column('min'); +#($pivot4row, $value4r, $unbounded4a) = $tableau1->find_pivot_row(3); +#($row4,$col4,$optimum4,$unbounded4) = $tableau1->find_next_pivot('min'); +#optimum is reached so can't go forward from here +#$basis4 = $tableau1->find_next_basis_from_pivot($row4, $col4); +#$tableau1->basis($basis4->value); +#$tableaudisplay3 = $tableau1->current_tableau; + +### experimental +# basis calculation +TEXT(" start with basis (6,7) and pivot on (1,1)"); + $basis = Set(6,7); + $basis = $basis - Set(6); + $basis = List( 1, $basis->value, ); +TEXT( "new basis is $basis"); + +############################################################## +# +# Text +# +# + +Context()->texStrings; +BEGIN_TEXT +matrix: \[ $m \quad $m2 \] $PAR +tableau: \[ $tableau_orig \] $PAR +original basis: $basis_orig $PAR +$HR + +\[\{ $tableau1->current_tableau\}\] +basis \{$tableau1->basis()\} ,\{List($tableau1->{basis})\} $PAR +basis matrix \[\{$tableau1->{M}-> + submatrix(rows=>[1..($tableau1->{m})],columns=>$tableau1->{basis})\}\]; +$PAR +new basis (2,3) +\[\{$tableau1->current_tableau(2,3)\}\] +basis \{$tableau1->basis()\}, \{join(",",@{$tableau1->{basis}})\}$PAR +basis matrix \[\{$tableau1->{M}->submatrix(rows=>[1..($tableau1->{m})],columns=>$tableau1->{basis})\}\]; +$PAR +new basis (3,2) +\[\{$tableau1->current_tableau(3,2)\}\] +basis \{$tableau1->basis()\}, \{join(",",@{$tableau1->{basis}})\}$PAR +basis matrix \[\{$tableau1->{M}->submatrix(rows=>[1..($tableau1->{m})],columns=>$tableau1->{basis})\}\]; +Should figure out why these don't reverse although internally they do. +$HR +submatrix (1,2) by (1,2) +\[\{$tableau1->{M}->submatrix(rows=>[1,2],columns=>[1,2])\}\]; +\[\[\{$tableau1->{B}\}\]; +submatrix (1,2) by (2,1) + +\[\{$tableau1->{M}->submatrix(rows=>[1,2],columns=>[2,1])\}\]; +\[\[\{$tableau1->{B}\}\]; +$HR +Figure out outputs: $PAR +\[\{$tableau1->current_tableau(1,7) \} \]$PAR +test0 ( \(\{join(" |", $tableau1->find_next_basis('min'))\}\) )$PAR +test1 ( \{join(" |", @{ List($tableau1->find_next_basis('min'))->data} )\} )$PAR +test2 ( \{ [$tableau1->find_next_basis('min')]->[0]->value \}, +\{ [$tableau1->find_next_basis('min')]->[1]->value \} )$PAR + +find next pivot column ( \{ List( $tableau1->find_pivot_column('min')) \} ) +----(\{ join ",", map {ref($_)} ( $tableau1->find_pivot_column('min')) \}) +$PAR +find next pivot row for col 2: ( \{ List($tableau1->find_pivot_row(2)) \} ) +----(\{ join ",", map {ref($_)} ( $tableau1->find_pivot_row(2)) \}) +$PAR +next pivot \( ( \{ List( $tableau1->find_next_pivot('min')) \} )\) +----(\{ join ",", map {ref($_)} ( $tableau1->find_next_pivot('min')) \}) +$PAR +next basis \( ( \{ List( $tableau1->find_next_basis('min')) \} ) \) +----(\{ join ",", map {ref($_)} ( $tableau1->find_next_basis('min')) \}) +$PAR +test pivot \( ( \{ List($tableau1->find_next_pivot('min')) \} ) \) +$PAR +END_TEXT + +Context()->normalStrings; + +ENDDOCUMENT(); + diff --git a/t/matrix_tableau_tests/tableau_test1.pg b/t/matrix_tableau_tests/tableau_test1.pg index cc0a503fe7..0983a57977 100644 --- a/t/matrix_tableau_tests/tableau_test1.pg +++ b/t/matrix_tableau_tests/tableau_test1.pg @@ -36,7 +36,7 @@ $m2 = $m->column_slice([3,2,1]); $list = $m->extract_rows_to_list(2,3); -$b = Matrix([1, 2, 3, 4]); +$b = Matrix([1, 2, 3, 4])->transpose; # make it an n by 1 matrix #TEXT($BR, "vector", $b->data->[1]); $c = Matrix([5, 6, 7]); $t = Tableau->new(A=>$m,b=>$b, c=>$c); @@ -49,11 +49,11 @@ $c_B2 = Value::Vector->new(map {$_->value} @$c_B); $c_4 = $t->current_tableau; -my $Badj = ($t->{B}->det) * ($t->{B}->inverse); +my $Badj = ($t->{current_basis_matrix}->det) * ($t->{current_basis_matrix}->inverse); my $current_tableau = $Badj * $t->{M}; # the A | S | obj | b $correction_coeff = ($c_B2*$current_tableau)->row(1); -$obj_row_normalized = ($t->{B}->det) *$t->{obj_row}; +$obj_row_normalized = ($t->{current_basis_matrix}->det) *$t->{obj_row}; $current_coeff = $obj_row_normalized-$correction_coeff ; TEXT("obj_row ", $t->{obj_row}, $BR ); @@ -71,7 +71,7 @@ original tableau is [`[$t->{M}]`] and basis is [$t->basis] -B is [`[$t->{B}]`] with determinant [$t->{B}->det] +B is [`[$t->{current_basis_matrix}]`] with determinant [$t->{current_basis_matrix}->det] the objective row is [@ $t->{obj_row} @] diff --git a/t/matrix_tableau_tests/tableau_test3.pg b/t/matrix_tableau_tests/tableau_test3.pg new file mode 100644 index 0000000000..705ad114e7 --- /dev/null +++ b/t/matrix_tableau_tests/tableau_test3.pg @@ -0,0 +1,843 @@ +DOCUMENT(); +loadMacros( +"PGstandard.pl", +"MathObjects.pl", +"parserPopUp.pl", +"unionLists.pl", +"MatrixReduce.pl", +#"AppletObjects.pl", +"PGessaymacros.pl", +"PGmatrixmacros.pl", +"LinearProgramming.pl", +"parserLinearInequality.pl", +"quickMatrixEntry.pl", +#"scaffold.pl", +"tableau.pl", +#"gage_matrix_ops.pl", +"PGinfo.pl", +"source.pl", +"PGcourse.pl", +); + +TEXT(beginproblem()); +TEXT($BEGIN_ONE_COLUMN); +$showPartialCorrectAnswers = 1; + +INITIALIZE_QUICK_MATRIX_ENTRY(); + + +############################################################## +# problem data +############################################################## + +# this can be changed (slightly at least) without affecting the behavior of the problem +# The choice of pivots is not yet automatic However + +# Your resources: +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + +#Hack to prevent domain conflict in answer. +Context()->variables->add(p1=>'Real',p2=>'Real'); +$objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); +# why can't the formula be defined within context "linearInequality"? + +Context("LinearInequality"); +Context()->variables->add(p1=>'Real',p2=>'Real'); +Context()->strings->add("Essay Answer" =>{}); +Context()->strings->add('Minimize'=>{},'Maximize'=>{}, "?"=>{}); +Context()->strings->add('Yes'=>{},'No'=>{}); +Context()->flags->set({tolType=>"absolute",tolerance=>.001}); +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +# $objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); + +our $context=Context(); + +############################################################## +# +# Setup +# +# + + + +$original_matrix = Matrix([ + [$bill_money_commitment, $steve_money_commitment, 1, 0, 0,0,0, $money_total], + [$bill_time_commitment,$steve_time_commitment, 0, 1, 0,0,0, $time_total], + [1,0,0,0,1,0,0,1], + [0,1,0,0,0,1,0,1], + [-$bill_profit, -$steve_profit, 0, 0, 0,0,1, 0] + ] +); +$toplabels = [qw(p1 p2 x3 x4 x5 x6 P b)]; +$sidelabels = [' ', qw(\text{cash} \text{hours} \text{p_1_bound} \text{p_2_bound} \text{obj_func}) ]; +$matrix1 = $original_matrix; + +############################################################################## +# utility subroutine for checking your answers +############################################################################## + +# size of constraint matrix with slack variables, objective column and constant column +$row_size =4; # n-1 n is constraints + objective_row +$param_size=2; # original number of parameters in the problem +$col_size = $param_size+$row_size+2; # m paramvars+slackvars+objcol+constant_col + +sub update_tableau { + my $original_matrix = shift; + my $basis = shift; + # $basis is a Set object + # last row is numbered $basis_size+1; + my $basis_size=@{$basis->data}; + my $basis_matrix = $original_matrix->submatrix(rows=>[1..$basis_size],columns=>$basis->data); + + my $obj_coeff = - ($original_matrix->row($basis_size+1)); + #obj function coeff corresponding to basis variables/columns + $cB = Matrix($obj_coeff->column_slice($basis->data)); + # get A | S |b portion of tableau and multiply by det(B)B^-1 + my $current_matrix= ($basis_matrix->det)* + ($basis_matrix->inverse)* + $original_matrix->row_slice(1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + return $current_tableau; +} +# the full tableau has one more row -- the objective function +sub display_tableau_state { + my ($original_matrix,$tableau, $matrix, $basis,$pivot)= @_; + $basis = $basis->sort; + my $basis_size = @{$basis->data} ; + my $basis_matrix = matrix_from_submatrix($original_matrix, rows=>[1..$basis_size],columns=>$basis->data); + my $normalized_tableau = $matrix; + my $reduced_matrix = matrix_from_matrix_rows($matrix,1..$basis_size); + warn "basis matrix is $basis_matrix, $reduced_matrix"; + my $normalized_reduced_matrix = ($basis_matrix->det)*$reduced_matrix; + # calculate current tableau by multiplying by the inverse of the basis + my $obj_coeff = - $original_matrix->row($basis_size+1); + my $cB = Matrix($obj_coeff->column_slice($basis->data)); + my $current_matrix = ($basis_matrix->det)*($basis_matrix->inverse)*matrix_from_matrix_rows($original_matrix,1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + $current_tableau = update_tableau($original_matrix, $basis); + + #DEBUG_MESSAGE("get state variables for $matrix using basis $basis"); + @statevars1 = get_tableau_variable_values($matrix, $basis); + # get z value + $statevars1 = ~~@statevars1; + $state = Matrix([[@statevars1]]); + +# return " +# +# pivot: \($pivot\) basis: \( $basis\) basis matrix \(". +# display_matrix_mm($basis_matrix)."\) $PAR +# tableau and normalized current tableau $PAR +# \(" . lp_display_mm($tableau) . "\)\(".lp_display_mm($normalized_tableau). " \) $PAR +# original matrix and current matrix calculated by Binverse*original $PAR +# \(". +# display_matrix_mm(matrix_from_matrix_rows($original_matrix,1..$basis_size)). +# "\) \(". +# display_matrix_mm($current_matrix). +# "\) $PAR +# objective coefficients = \( $obj_coeff \) $PAR +# basis objective_coefficients = \($cB\) $BR +# current matrix rows \(". +# join(" ", @{$current_matrix->extract_rows}). +# "\)$PAR +# new objective row cB * B^-1 * M -c= \($new_obj_row\)$PAR +# new tableau \[". +# lp_display_mm($current_tableau). +# "\] $PAR +# state: \(" . join(", ", @$state) . "\) $PAR +# " + return " + + pivot: \($pivot\) basis: \( $basis\) $PAR basis matrix \(". + display_matrix_mm($basis_matrix)."\) $PAR + tableau and normalized current tableau $PAR + \(" . lp_display_mm($tableau) . "\)\(".lp_display_mm($normalized_tableau). " \) $PAR + original matrix and current matrix calculated by Binverse*original $PAR + \(". + display_matrix_mm(matrix_from_matrix_rows($original_matrix,1..$basis_size)). + "\) \(". + display_matrix_mm($current_matrix). + "\) $PAR + objective coefficients = \( $obj_coeff \) $PAR + basis objective_coefficients = \($cB\) $PAR + new tableau \[". + lp_display_mm($current_tableau). + "\] $PAR + state: \($state)\) $PAR + " +} + +# here is the outline of the worksheet project +# BEGIN_TEXT +# \{BeginList("OL")\} +# $ITEM Set up problem +# $ITEM Write tableau +# $ITEM Solve using simplex method +# $ITEM Create dual problem +# $ITEM Write tableau for dual problem +# $ITEM Set up LOP for auxiliary method to find feasible solution +# $ITEM Solve the auxiliary problem to find a feasible basic solution. +# $ITEM Finish solving the dual problem using simplex method +# $ITEM Compare answers to the primary problem and the dual problem. +# \{EndList("OL")\} +# +# END_TEXT + + +######################### +# start sections +######################### + + + +######################################################################################################### +# Section 1 +# section 1 (1essay) +######################################################################################################### +Context($context); +$constraint1 = Compute("${bill_money_commitment}p1 + ${steve_money_commitment}p2 <= $money_total"); +$constraint2 = Compute("${bill_time_commitment}p1 + ${steve_time_commitment}p2 <= $time_total"); +$constraint3 = Compute("p1<=1"); +$constraint4 = Compute("p2<=1"); +$constraints = List( $constraint1, $constraint2, $constraint3, $constraint4 ); +$popup1 = PopUp([qw(? Maximize Minimize)],'Maximize'); +Context($context); +# defined above? wtf? # $objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); +Context()->texStrings; + +BEGIN_TEXT +$PAR This is the same as an earlier problem but the numbers have +been changed slightly. How does this change the result? +$HR $PAR 1.$PAR +You have $DOLLAR$money_total to invest. $PAR +Two of your friends, Bill and Steve, have offered you an opportunity to become +a partner in two different entrepreneurial ventures, one planned by each friend. In both cases, this +investment would involve expending some of your time next summer as well as putting up cash. +$PAR Becoming +a full partner in Bill's venture will require an investment of $DOLLAR$bill_money_commitment + and $bill_time_commitment hours and your estimated profit +at the end of the summer (ignoring the value of your time) would be +$DOLLAR$bill_profit (plus your investment money back). +$PAR +The corresponding figures for Steve's venture are +$DOLLAR$steve_money_commitment and $steve_time_commitment hours with an estimated profit of $DOLLAR$steve_profit. +Both friends are flexible and would +allow you to come in at any fraction of a full partnership you would like. If you choose a fraction +of a full partnership, all the above figures given for a full partnership (money investment, +time investment, and your profit) would be multiplied by this same fraction. +$PAR +Because you were looking for an interesting summer job anyway (maximum of $time_total hours), you have decided to participate in +one or both friends' ventures in whichever combination would maximize your total estimated profit. +$PAR +Write a linear program to help you determine the correct fractions. Use p1 for the fraction of +a full partnership with Bill and p2 for the fraction of a full partnership with Steve. These +fractions should each be less than or equal to one but don't have to sum to one. If you had enough money and +time you could invest fully in both ventures. +$PAR + +\{ANS($popup1->cmp), $popup1->menu\} the profit + +\(P = \)\{Context($context),ANS($objfun1->cmp), ans_rule(40)\} + +$PAR subject to the following constraints. Separate the constraints by commas. + +\{ ANS($constraints->cmp->withPostFilter( + linebreak_at_commas() +)), ans_box(10,80) \} +$PAR +END_TEXT + + +BEGIN_TEXT + +SOLUTION $PAR + +\(p_1\) is the fraction to invest in Bill's venture. +\(p_2 = \) fraction of Steve's venture to invest in. +$PAR + +Maximize \(P = $objfun1\) subject to the constraints: + $PAR +\[ +\begin{aligned} +$constraint1 & \\ +$constraint2 & \\ +$constraint3 & \\ +$constraint4 & \\ +\end{aligned} +\] + +END_TEXT + +######################################################################################################### +# Section 2 -- set up tableau +######################################################################################################### +# get information on current state +$tableau1 = $matrix1->wwMatrix->array_ref; # translate to the array reference +$basis1 = Set(3,4,5,6); +@statevars1 = get_tableau_variable_values($matrix1, $basis1); +# get z value +$statevars1 = ~~@statevars1; +$state1 = Matrix([[@statevars1]]); + +$matrix1->{top_labels}=$toplabels; +Context()->texStrings; + +BEGIN_TEXT + +Write the matrix/tableau representing the linear optimization problem above. Use +the convention that the objective function is listed on the bottom row and the coefficient in +front of the profit \(P\) is \(1\) or equivalently in the form \( -ax_1 -bx_2 +z = 0 \) +$PAR +We'll use x3 for the slack variable for the money constraint, x4 for the time constraint +slack variable and x5 and x6 for the slack variables for the contraints on p1 and p2. +$PAR +\{MATRIX_ENTRY_BUTTON($matrix1)\} +$PAR +\{ANS($matrix1->cmp), $matrix1->ans_array()\} +$PAR +END_TEXT + +BEGIN_TEXT +SOLUTION:$PAR + +\[ \{lp_display_mm($matrix1, top_labels=>$toplabels).side_labels($sidelabels)\} \] + +$PAR +First display of state 1 +$PAR +\{display_tableau_state($original_matrix,$tableau1, $matrix1, $basis1, $pivot)\} +$PAR +END_TEXT + + +Context()->normalStrings; + +######################################################################################################### +# Section 3 -- pivots +######################################################################################################### + +# ######################################################################################################### +# # first pivot +# ######################################################################################################### +# +# Context($context); +$pivot1 = Point("(3,1)"); +# the new basis vectors are [1,3,4,6] <- [3,4,5,6] +($tableau2,$basis2, $statevars2) = lp_basis_pivot($tableau1,$basis1,$pivot1); + $B2 = matrix_extract_submatrix($original_matrix,rows=>[1..$row_size],columns=>$basis2->data); +# BEGIN_TEXT +# tableau2 = $tableau2 $BR +# basis2 = $basis2 $BR +# statevars2 = $statevars2 $BR +# original_matrix = $original_matrix $BR +# B2 = $B2 $BR +# row data:\{ join(" ", @{[1..$row_size]})\} $BR +# basis data (columns: \{join(" ", @{ $basis2->data })\} $BR +# END_TEXT +$matrix2 = ($B2->det)*Matrix($tableau2); +$matrix2->{top_labels}=$toplabels; +$state2=Matrix([[@$statevars2]]); + + +BEGIN_TEXT +SOLUTION:$PAR + +\[ \{lp_display_mm($matrix2, top_labels=>($matrix2->{top_labels})).side_labels($sidelabels)\} \] + +$PAR +Display of state 2 $PAR +\{display_tableau_state($original_matrix,$tableau2, $matrix2, $basis2, $pivot1)\} + +END_TEXT + +# BEGIN_TEXT +# +# Some experiments +# extract rows +# \{ @temp = @{$matrix2->extract_rows} \} +# \{pop @temp \} +# $PAR +# first element \{ref($temp[0])\} +# $PAR +# \{ join(" ", @temp, Matrix([1,2,3,4,5,6,7,8])) \}$PAR +# +# \{$temp2 = Matrix(@temp, Matrix([1,2,3,4,5,6,7,8])) \} +# $PAR +# new matrix is \[ \{lp_display_mm($temp2, top_labels=>$toplabels).side_labels($sidelabels)\} \] +# +# +# END_TEXT +# ######################################################################################################### +# # second pivot +# ######################################################################################################### +$pivot2 = Point("(1,2)"); +# the new basis vectors are [1,2,4,6] <- [1,3,4,6] +($tableau3,$basis3, $statevars3) = lp_basis_pivot($tableau2,$basis2,$pivot2); +$B3 = matrix_from_submatrix($original_matrix,rows=>[1..$row_size],columns=>$basis3->data); +$matrix3 = ($B3->det)*Matrix($tableau3); +$matrix3->{top_labels}=$toplabels; +$state3=Matrix([[@$statevars3]]); + +######################################################################################################### +# third pivot +######################################################################################################### +$pivot3 = Point("(2,5)"); +# the new basis vectors are [1,2,5,6] <- [1,2,4,6] +($tableau4,$basis4, $statevars4) = lp_basis_pivot($tableau3,$basis3,$pivot3); +$B4 = matrix_from_submatrix($original_matrix,rows=>[1..$row_size],columns=>$basis4->data); +$matrix4 = ($B4->det)*Matrix($tableau4); +$matrix4->{top_labels}=$toplabels; +$state4=Matrix([[@$statevars4]]); + +$p1 = $state4->element(1,1); +$p2 = $state4->element(1,2); +$profit = $state4->element(1,7); + +######################################################################################################### +# debugging display code +######################################################################################################### +BEGIN_TEXT +pivot1 $pivot1 $PAR +\{display_tableau_state($original_matrix,$tableau2, $matrix2, $basis2, $pivot1)\} +pivot2 $pivot2 $PAR +\{display_tableau_state($original_matrix,$tableau3, $matrix3, $basis3, $pivot2)\} +pivot3 $pivot3 $PAR +\{display_tableau_state($original_matrix,$tableau4, $matrix4, $basis4, $pivot3)\} + +$PAR +end of three pivots +END_TEXT +# # +# # +# # ######################################################################################################### +# # # section 3 answers (3,4,5,6,7,8, 9, 10) +# # ######################################################################################################### +Context()->texStrings; + +BEGIN_TEXT +$PAR +What is the initial state? $PAR + \((p_1,p_2,x_3,x_4,x_5,x_6,P)=\): \{ANS($state1->cmp),$state1->ans_array(4)\} $PAR + +Using the convention that one removes the "first" nonbasic column that will increase profits +(the convention that Hurlbert uses) find the first pivot location. (e.g. (2,3) for line 2, column 3). +and perform the pivot operation. (You can use +\{htmlLink("http://people.vcu.edu/~ghurlbert/websim/sim.html","websim")\} to do the calculation. ); +$PAR + +Pivot entry: use parentheses as in (row, column): \{ANS($pivot1->cmp),$pivot1->ans_rule \}. +$PAR +The matrix has been multiplied by a constant so that all values are integers. This is +the same result as you will obtain using websim. +$PAR +\{MATRIX_ENTRY_BUTTON($matrix2)\} +$PAR +$PAR +\{ANS($matrix2->cmp()), $matrix2->ans_array(6)\} + +$PAR +Current vector (values): \{ANS($state2->cmp),$state2->ans_array(6)\} $PAR +$HR +Next pivot entry: \{ANS($pivot2->cmp),$pivot2->ans_rule \}. + +$PAR +\{MATRIX_ENTRY_BUTTON($matrix3)\} +$PAR +\{ANS($matrix3->cmp()), $matrix3->ans_array(6)\} +$HR +Next pivot entry: \{ANS($pivot3->cmp),$pivot3->ans_rule \} +results in the tableau: +$PAR +\{MATRIX_ENTRY_BUTTON($matrix4)\} +$PAR $PAR +\{ANS($matrix4->cmp()), $matrix4->ans_array(6)\} +$HR +At this point we are at a local (and therefore a global) maximum and pivoting +will not increase the profit. +$PAR +Your final state is \{ANS($state4->cmp),$state4->ans_array(6)\} $PAR +(Make sure you are getting all the zeros in the answer! The numbers are quite large +and overflow the blanks. -- you can make the table blanks wider +in websim by right clicking on upper left corner of a table. ) + +$PAR What percentage of Bill's venture do you invest in? +\{ANS($p1->cmp), $p1->ans_rule\} +$PAR What percentage of Steve's venture do you invest in? +\{ANS($p2->cmp), $p2->ans_rule\} +$PAR What is your expected profit? +$DOLLAR\{ANS($profit->cmp),$profit->ans_rule\}. + + +$PAR + +END_TEXT + +BEGIN_SOLUTION +Before pivoting the state is \($state1\) $PAR +Your first pivot should be on \($pivot1\) since that is the left most +column that +will increase the profit and the first row has the most limiting ratio. +$PAR +\( \{lp_display_mm( [$matrix2->value],top_labels=>$toplabels )\} \) $PAR +$PAR + +$PAR + \( \{lp_display_mm([$matrix3->value], top_labels=>$toplabels )\}\) +\( \{side_labels( qw(\text{cash} \text{hours} \text{profits} ) ) \} \) +$PAR + + +one more + \( \{lp_display_mm([$matrix4->value],top_labels=>$toplabels )\} \) + $PAR +At this point you are done since changing any of the non basic variables will only +decrease the profit. +END_SOLUTION +Context()->normalStrings; +# + +# ######################################################################################################### +# # Section 4 dual problem +# ######################################################################################################### + + +Context()->variables->add(y1=>'Real',y2=>'Real', y3=>'Real', y4=>'Real', w=>'Real',y0=>'Real'); +$dual_constraint1 = Compute("${bill_money_commitment}y1 +${bill_time_commitment}y2 +y3 >=$bill_profit"); +$dual_constraint2 = Compute("${steve_money_commitment}y1 +${steve_time_commitment}y2 +y4 >=$steve_profit"); +$dual_constraints= List($dual_constraint1, $dual_constraint2); + +$dual_objfun = Formula("${money_total}y1 +${time_total}y2 + y3 +y4"); +$popupmaxmin = PopUp(["?","Maximize", "Minimize"], "Minimize"); + +Context()->texStrings; + +BEGIN_TEXT +Construct the dual problem for the linear optimization problem above. The first goal is to calculate +an upper bound for the possible profit in the LOP using linear combinations of the inequality constraints. +Then formulate the search for the best of these upper bounds in such a way that it +becomes a new LOP -- the dual problem. Use variables \(y_1,y_2,y_3,y_4\) to create linear combinations +of the constraints on money, time and the total probabilities \(p_1\) and \(p_2\) respectively, +and create a linear function \(w = Ay_1 + By_2 +Cy_3 +Dy_4\) which guarantees that \(w\) will +be larger than any profit one could make. What constraints must \(y_1\dots y_4\) satisfy? +$PAR +\{ANS($dual_constraints->cmp->withPostFilter( + linebreak_at_commas() +)),ans_box(4,80)\} +$PAR The objective function would be: + \(w = \) \{ANS($dual_objfun->cmp), $dual_objfun->ans_rule(50)\} +$PAR +To get the best possible estimate for the profit in the original problem +you would want to \{ANS($popupmaxmin->cmp),$popupmaxmin->menu\} \(w\) +Why? +$PAR +\{ANS(essay_cmp()), essay_box(5,80)\} +END_TEXT + +BEGIN_SOLUTION +The profit is less than the minimum of \( w = $dual_objfun \) subject to + +\[\begin{aligned} +$dual_constraint1 & \\ +$dual_constraint2 & +\end{aligned} +\] +The \(y_i\) are the coefficients used to add up the constraint inequalities of the primary problem +so as to estimate an upper bound for the profit. If the \(y_i\) are positive values satisfying +the constraints then \(w\) will be greater than or equal to the profit. We get the most precise +estimate by finding values of \(y_i\), satisfying the constraints, which give the smallest +possible value for \( w\). + +END_SOLUTION +Context()->normalStrings; + + +# ######################################################################################################### +# # Section 5 set up dual constraints +# ######################################################################################################### + +Context($context); +$dualtableau1 = [[-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1, 0, 0,-$bill_profit], + [-$steve_money_commitment, -$steve_time_commitment,0,-1, 0, 1, 0,-$steve_profit], + [${money_total},${time_total},1, 1, 0,0, 1,0]]; +$dualmatrix1 = Matrix($dualtableau1); +$dualtoplabels = [qw(y1 y2 y3 y4 y5 y6 v b)]; +$dualmatrix1->{top_labels}=$dualtoplabels; +$dualtableau1_string = lp_display_mm($dualtableau1,top_labels=>$dualtoplabels); + +$popup = PopUp(["?","Yes", "No"], "No"); + +$dual_row_size =2; # figure these out automatically +$dual_col_size =$col_size; +BEGIN_TEXT +dualmatrix1 size \{join(",", $dualmatrix1->dimensions)\} $BR +$dual_row_size, $dual_col_size +END_TEXT +Context()->texStrings; + +BEGIN_TEXT +$PAR +For this first effort let's rewrite the equations so that they are in "standard" form, meaning +that the inequalities are all less than or equal to and that the goal is to maximize \(v = -w\). +This will probably mean that you will have to rewrite the natural way in which you set up +the problem above. Some coefficients will change sign. +We will still be using the convention that +the coefficient of \(v\) will be \(1\). + +$PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix1)\} +\{ANS($dualmatrix1->cmp()), $dualmatrix1->ans_array \} +$PAR +(This tableau is for maximizing \(v\), aka \(-w\) -- there are two minus sign switches. ) +$PAR +Is there a natural feasible solution to this problem? In other words does setting the +problem parameters equal to 0 provide a feasible solution? +\{ANS($popup->cmp),$popup->menu \} +$PAR +Explain why or why not and what are your options for getting started. +\{ ANS(essay_cmp()), essay_box(10,80) \} +$PAR +END_TEXT + +BEGIN_SOLUTION +\[ $dualtableau1_string \] +$PAR +Setting the parameters \(y_1, y_2\) to zero is not a feasible solution. +Options for finding a first feasible solution include guessing (not a bad choice +for a problem this small), using prior knowledge (e.g. a known optimal solution to +a similar problem), creating an auxiliary problem to find a feasible point, +and the "shortcut method" which is an accelerated version of the "auxiliary method". +END_SOLUTION + + +# ######################################################################################################### +# # Begin section 6 -- solve the tableau using simplex method +# # section 6 (13matrix, 14essay) +# # First we have to find a feasible solution -- phase 1 +# # We'll use the auxiliary method -- adding an extra slack variable. +# ######################################################################################################### + + +Context($context); +$dualtableau2aux = [[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1,0,0,0,-$bill_profit], + [-1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, 0,1,0,0,-$steve_profit], + [1,0,0,0,0,0,0,0,1,0]]; + +# [0,$money_total,$time_total,1,1,0,0,1,0,0]]; +$dualmatrix2aux = Matrix($dualtableau2aux); +$original_dual_aux_matrix=$dualmatrix2aux; +$dualtoplabelsphase1 = [qw(y0 y1 y2 y3 y4 y5 y6 v z b)]; +$dualmatrix2aux->{top_labels}= $dualtoplabelsphase1; +$dual_constraint1phase1 = Compute("-${bill_money_commitment}y1 -${bill_time_commitment}y2 -y3 <=-$bill_profit+y0"); +$dual_constraint2phase1 = Compute("-${steve_money_commitment}y1 -${steve_time_commitment}y2 -y4 <=-$steve_profit+y0"); +$dualconstraintsphase1 = List($dual_constraint1phase1,$dual_constraint2phase1); +$dualpivot2aux=Point(1,1); +$dualbasis2aux=Set(6,7); +######################################################################################################### +# section 6 answers (15essay 16matrix) +######################################################################################################### +Context()->texStrings; + +BEGIN_TEXT +Using the auxiliary method write the LOP for finding +the first feasible point. Write +the new constraints to be used by auxiliary method. $PAR +\{ ANS($dualconstraintsphase1->cmp),ans_box(3,80) \} +$PAR +We will continue to use the convention that we are maximizing +the auxiliary function so there will be some minus signs that need to be taken into account. +We'll let \(z = -y_0\) and try to maximize \(z\). +$PAR +Now construct the tableau you'll use for the auxiliary method. +Add a first column for the extra slack variable and a next to the last row +for the new objective function value \(z\) and write the new tableau. +The next to the last row will hold the objective function for \(z\). +The last row will hold the original +dual objective function which we'll just carry along. +$PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix2)\} +\{ANS($dualmatrix2aux->cmp()), $dualmatrix2aux->ans_array \} +$PAR +END_TEXT +######################################################################################################### +# debugging display code +######################################################################################################### +BEGIN_TEXT + +Display dualtableau2aux +\($dualtableau2aux\), matrix \($dualmatrix2aux\) +dualbasis \($dualbasis2aux\) and \($dualpivot2aux\) +$PAR +\{display_tableau_state($original_dual_aux_matrix,$dualtableau2aux, $dualmatrix2aux, $dualbasis2aux, $dualpivot2aux)\} +$PAR +END_TEXT + +Context()->normalStrings; + +Context()->texStrings; +BEGIN_SOLUTION +New objective: Maximize \(z = -y_0 \) subject to +\[\begin{align} +$dual_constraint1phase1&\\ +$dual_constraint2phase1 &\\ +\end{align} +\] with all variables non-negative. +$PAR +If \(y_0\) is large enough this always has a feasible solution with +\(y_1=y_2=0\). + +\[ \{lp_display_mm($dualtableau2aux) \} \] +$PAR +END_SOLUTION +Context()->normalStrings; + + + +# ######################################################################################################### +# # Section 7 final comparisons of primary and dual results. +# ######################################################################################################### + ($dualtableau3aux,$dualbasis3aux, $dualstatevars3aux) = lp_basis_pivot($dualtableau2aux,$dualbasis2aux,$dualpivot2aux); + $B3 = matrix_from_submatrix($original_dual_aux_matrix,rows=>[1..$dual_row_size],columns=>($dualbasis3aux->data)); +$dualmatrix3aux = ($B3->det)*Matrix($dualtableau3aux); + $dualmatrix3aux->{top_labels}=$dualtoplabelsphase1; + $dualstate3aux=Matrix([[@$dualstatevars3aux]]); +######################################################################################################### +# debugging display code +######################################################################################################### +BEGIN_TEXT +$PAR +Debugging display code dualtableau3aux. This is the matrix +after pivoting on (1,1). $PAR + +\{display_tableau_state($original_dual_aux_matrix,$dualtableau3aux, $dualmatrix3aux, $dualbasis3aux, $dualpivot3aux)\} +$PAR +END_TEXT +# +# $pivot4 = Point("(1,1)"); +# $dualtableau3 = lp_pivot($dualtableau2,0,0); +# $dualmatrix3 = Matrix($dualtableau3); +# $dualmatrix3->{top_labels}=$dualtoplabelsphase1; +# $dualtableau3_string = lp_display_mm($dualtableau3); +# $z_initial=-$bill_profit; # a hint +# # this is the initial value of the z=-y0 variable we are maximizing? +# +# # now start to maximize z +# $pivot5 = Point("(2,2)"); +# $dualtableau4 = lp_pivot($dualtableau3,1,1); +# $denom1 = 2000; +# $dualmatrix4 = $denom1*Matrix($dualtableau4); +# $dualmatrix4->{top_labels}=$dualtoplabelsphase1; +# $dualtableau4_string = lp_display_mm($dualmatrix4); +# $pivot6 = Point("(1,3)"); +# $dualtableau5 = lp_pivot($dualtableau4,0,2); +# +# ## FIXME entries close to zero in a matrix are not being compared properly. +# ## FIXME tolerance and tolType are being ignored +# +# #$dualtableau5->[2][3]=0; +# #$dualtableau5->[2][5]=0; +# +# $denom2 = 1.3E6; +# $dualmatrix5 = $denom2*Matrix($dualtableau5); +# $dualmatrix5->{top_labels}=$dualtoplabelsphase1; +# $dualtableau5_string = lp_display_mm($dualmatrix5); +# $state5 = Matrix([[ 6.4615, 0.423, 0, 0, -6415.38]]); +# $popup5 = PopUp(["?","Yes", "No"], "No"); +# ######################################################################################################### +# # Section 7 answers 17,18matrix,19,20matrix,21,22matrix,23, 24,25essay +# ######################################################################################################### +# +# +# +# Context()->texStrings; +# +# BEGIN_TEXT +# The first pivot, which will make the right hand side entries positive is +# \{ANS($pivot4->cmp),$pivot4->ans_rule \}. $PAR +# Recall that in the case of tie we take entry with the least index +# (i.e. left most or upper most. ) +# $PAR +# +# +# The resulting tableau is $PAR +# \{MATRIX_ENTRY_BUTTON($dualmatrix3)\} +# \{ANS($dualmatrix3->cmp), $dualmatrix3->ans_array() \} +# +# The value of \(z=-y_0\) is \($z_initial\). +# $PAR +# The next pivot, following the simplex method to maximize \(z\), is +# \{ANS($pivot5->cmp), $pivot5->ans_rule\}. +# $PAR +# (Many of the following answers have lots of zeros. You can use the shortcut +# 1E3 to stand for \(1\times 10^3\). ) +# $PAR +# Notice that because of the zero on the right hand side none of the state variables change. We had three +# hyperplanes intersecting at a point and we have changed our mind about which +# of those three we consider basic. The new tableau is +# $PAR +# \{MATRIX_ENTRY_BUTTON($dualmatrix4)\} +# \{ANS($dualmatrix4->cmp), $dualmatrix4->ans_array() \} +# $PAR +# The next pivot \{ ANS($pivot6->cmp), $pivot6->ans_rule \} leads to $PAR +# \{MATRIX_ENTRY_BUTTON($dualmatrix5)\} +# \{#FIXME -- these tolerance leveals are being ignored \} +# \{ANS($dualmatrix5->cmp(zeroLevel=>.1, zeroLevelTol=>.1)), $dualmatrix5->ans_array() \} +# +# $PAR and we notice that now \(z=-y_0=0\) so we have found a basic feasible solution to our +# original dual problem. The variables \(y_1,y_2,y_3,y_4,v\) for this +# solution are +# $PAR +# +# \{ANS($state5->cmp),$state5->ans_array\} +# $PAR +# Do we need to continue to optimize the value for \(v\)? +# \{ANS($popup5->cmp), $popup5->menu()\} Why? $PAR +# \{ANS(essay_cmp), essay_box(3,80)\} +# $PAR +# Compare this answer \(v^*= -w^*\) to the dual problem to the optimal value \(P^*\)for the primary problem. +# The problems' goals were to maximize \(P\) and to minimize \(w\). +# $PAR +# +# END_TEXT +# +# BEGIN_SOLUTION +# The first pivot is \($pivot4\) which makes the right hand side entries positive. +# $PAR +# \[$dualtableau3_string\] +# $PAR +# The next pivot is \($pivot5\). It follows the simplex rule of choosing the row +# with the most restrictive ratio -- in this case zero. The result is simply to +# choose a new representation of the same point -- no change in state takes place. +# $PAR\[$dualtableau4_string\] +# $PAR +# The final pivot (as it turns out) is \($pivot6\). At this point we are done with +# phase 1 because the value of \(z=-y_0=0\) so we have maximized \(z\) and +# minimized \(y_0\) to \(0\) +# $PAR\[$dualtableau5_string\] +# $PAR +# At this point we notice that the coefficients in the last row are such that +# we cannot increase the value of \(v\) any further. \(v\) is at its maximum +# and \(w=-v\) is at its minimum. The minimum value of \(w=6000\) which is the +# same as the maximum profit that we found in the first example. +# END_SOLUTION +# Context()->normalStrings; +# Section::End(); +# TEXT($END_ONE_COLUMN); +# Scaffold::End(); +# } # END SECTION 7 +ENDDOCUMENT(); + diff --git a/t/matrix_tableau_tests/tableau_test4.pg b/t/matrix_tableau_tests/tableau_test4.pg new file mode 100644 index 0000000000..02448d246f --- /dev/null +++ b/t/matrix_tableau_tests/tableau_test4.pg @@ -0,0 +1,512 @@ +DOCUMENT(); +loadMacros( +"PGstandard.pl", +"MathObjects.pl", +"parserPopUp.pl", +"unionLists.pl", +"MatrixReduce.pl", +#"AppletObjects.pl", +"PGessaymacros.pl", +"PGmatrixmacros.pl", +"LinearProgramming.pl", +"parserLinearInequality.pl", +"quickMatrixEntry.pl", +#"scaffold.pl", +"tableau.pl", +#"gage_matrix_ops.pl", +"PGinfo.pl", +"source.pl", +"PGcourse.pl", +); + +TEXT(beginproblem()); +TEXT($BEGIN_ONE_COLUMN); +$showPartialCorrectAnswers = 1; + +INITIALIZE_QUICK_MATRIX_ENTRY(); + + +############################################################## +# problem data +############################################################## + +# this can be changed (slightly at least) without affecting the behavior of the problem +# The choice of pivots is not yet automatic However + +# Your resources: +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + +#Hack to prevent domain conflict in answer. +Context()->variables->add(p1=>'Real',p2=>'Real'); +$objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); +# why can't the formula be defined within context "linearInequality"? + +Context("LinearInequality"); +Context()->variables->add(p1=>'Real',p2=>'Real'); +Context()->strings->add("Essay Answer" =>{}); +Context()->strings->add('Minimize'=>{},'Maximize'=>{}, "?"=>{}); +Context()->strings->add('Yes'=>{},'No'=>{}); +Context()->flags->set({tolType=>"absolute",tolerance=>.001}); +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +# $objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); + +our $context=Context(); + +############################################################## +# +# Setup +# +# + + + +############################################################################## +# utility subroutine for checking your answers +############################################################################## + +sub update_tableau { + my $original_matrix = shift; + my $basis = shift; + # $basis is a Set object + # last row is numbered $basis_size+1; + my $basis_size=@{$basis->data}; + my $basis_matrix = $original_matrix->submatrix(rows=>[1..$basis_size],columns=>$basis->data); + + my $obj_coeff = - ($original_matrix->row($basis_size+1)); + #obj function coeff corresponding to basis variables/columns + $cB = Matrix($obj_coeff->column_slice($basis->data)); + # get A | S |b portion of tableau and multiply by det(B)B^-1 + my $current_matrix= ($basis_matrix->det)* + ($basis_matrix->inverse)* + $original_matrix->row_slice(1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + return $current_tableau; +} +# the full tableau has one more row -- the objective function +sub display_tableau_state { + my ($original_matrix,$tableau, $matrix, $basis,$pivot)= @_; + $basis = $basis->sort; + my $basis_size = @{$basis->data} ; + my $basis_matrix = matrix_from_submatrix($original_matrix, rows=>[1..$basis_size],columns=>$basis->data); + my $normalized_tableau = $matrix; + my $reduced_matrix = matrix_from_matrix_rows($matrix,1..$basis_size); + warn "basis matrix is $basis_matrix, $reduced_matrix"; + my $normalized_reduced_matrix = ($basis_matrix->det)*$reduced_matrix; + # calculate current tableau by multiplying by the inverse of the basis + my $obj_coeff = - $original_matrix->row($basis_size+1); + my $cB = Matrix($obj_coeff->column_slice($basis->data)); + my $current_matrix = ($basis_matrix->det)*($basis_matrix->inverse)*matrix_from_matrix_rows($original_matrix,1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + $current_tableau = update_tableau($original_matrix, $basis); + + #DEBUG_MESSAGE("get state variables for $matrix using basis $basis"); + @statevars1 = get_tableau_variable_values($matrix, $basis); + # get z value + $statevars1 = ~~@statevars1; + $state = Matrix([[@statevars1]]); + + + return " + + pivot: \($pivot\) basis: \( $basis\) $PAR basis matrix \(". + display_matrix_mm($basis_matrix)."\) $PAR + tableau and normalized current tableau $PAR + \(" . lp_display_mm($tableau) . "\)\(".lp_display_mm($normalized_tableau). " \) $PAR + original matrix and current matrix calculated by Binverse*original $PAR + \(". + display_matrix_mm(matrix_from_matrix_rows($original_matrix,1..$basis_size)). + "\) \(". + display_matrix_mm($current_matrix). + "\) $PAR + objective coefficients = \( $obj_coeff \) $PAR + basis objective_coefficients = \($cB\) $PAR + new tableau \[". + lp_display_mm($current_tableau). + "\] $PAR + state: \($state)\) $PAR + " +} + + +# ######################################################################################################### +# # Section 4 dual problem +# ######################################################################################################### + + +Context()->variables->add(y1=>'Real',y2=>'Real', y3=>'Real', y4=>'Real', w=>'Real',y0=>'Real'); +$dual_constraint1 = Compute("${bill_money_commitment}y1 +${bill_time_commitment}y2 +y3 >=$bill_profit"); +$dual_constraint2 = Compute("${steve_money_commitment}y1 +${steve_time_commitment}y2 +y4 >=$steve_profit"); +$dual_constraints= List($dual_constraint1, $dual_constraint2); + +$dual_objfun = Formula("${money_total}y1 +${time_total}y2 + y3 +y4"); +$popupmaxmin = PopUp(["?","Maximize", "Minimize"], "Minimize"); + +Context()->texStrings; + +BEGIN_TEXT +Construct the dual problem for the linear optimization problem above. The first goal is to calculate +an upper bound for the possible profit in the LOP using linear combinations of the inequality constraints. +Then formulate the search for the best of these upper bounds in such a way that it +becomes a new LOP -- the dual problem. Use variables \(y_1,y_2,y_3,y_4\) to create linear combinations +of the constraints on money, time and the total probabilities \(p_1\) and \(p_2\) respectively, +and create a linear function \(w = Ay_1 + By_2 +Cy_3 +Dy_4\) which guarantees that \(w\) will +be larger than any profit one could make. What constraints must \(y_1\dots y_4\) satisfy? +$PAR +\{ANS($dual_constraints->cmp->withPostFilter( + linebreak_at_commas() +)),ans_box(4,80)\} +$PAR The objective function would be: + \(w = \) \{ANS($dual_objfun->cmp), $dual_objfun->ans_rule(50)\} +$PAR +To get the best possible estimate for the profit in the original problem +you would want to \{ANS($popupmaxmin->cmp),$popupmaxmin->menu\} \(w\) +Why? +$PAR +\{ANS(essay_cmp()), essay_box(5,80)\} +END_TEXT + +BEGIN_SOLUTION +The profit is less than the minimum of \( w = $dual_objfun \) subject to + +\[\begin{aligned} +$dual_constraint1 & \\ +$dual_constraint2 & +\end{aligned} +\] +The \(y_i\) are the coefficients used to add up the constraint inequalities of the primary problem +so as to estimate an upper bound for the profit. If the \(y_i\) are positive values satisfying +the constraints then \(w\) will be greater than or equal to the profit. We get the most precise +estimate by finding values of \(y_i\), satisfying the constraints, which give the smallest +possible value for \( w\). + +END_SOLUTION +Context()->normalStrings; + + +# ######################################################################################################### +# # Section 5 set up dual constraints +# ######################################################################################################### + +Context($context); +$dualtableau1 = [[-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1, 0, 0,-$bill_profit], + [-$steve_money_commitment, -$steve_time_commitment,0,-1, 0, 1, 0,-$steve_profit], + [${money_total},${time_total},1, 1, 0,0, 1,0]]; +$dualmatrix1 = Matrix($dualtableau1); +$dualtoplabels = [qw(y1 y2 y3 y4 y5 y6 v b)]; +$dualmatrix1->{top_labels}=$dualtoplabels; +$dualtableau1_string = lp_display_mm($dualtableau1,top_labels=>$dualtoplabels); + +$popup = PopUp(["?","Yes", "No"], "No"); + +$dual_row_size =2; # figure these out automatically +$dual_col_size =$col_size; +BEGIN_TEXT +dualmatrix1 size \{join(",", $dualmatrix1->dimensions)\} $BR +$dual_row_size, $dual_col_size +END_TEXT +Context()->texStrings; + +BEGIN_TEXT +$PAR +For this first effort let's rewrite the equations so that they are in "standard" form, meaning +that the inequalities are all less than or equal to and that the goal is to maximize \(v = -w\). +This will probably mean that you will have to rewrite the natural way in which you set up +the problem above. Some coefficients will change sign. +We will still be using the convention that +the coefficient of \(v\) will be \(1\). + +$PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix1)\} +\{ANS($dualmatrix1->cmp()), $dualmatrix1->ans_array \} +$PAR +(This tableau is for maximizing \(v\), aka \(-w\) -- there are two minus sign switches. ) +$PAR +Is there a natural feasible solution to this problem? In other words does setting the +problem parameters equal to 0 provide a feasible solution? +\{ANS($popup->cmp),$popup->menu \} +$PAR +Explain why or why not and what are your options for getting started. +\{ ANS(essay_cmp()), essay_box(10,80) \} +$PAR +END_TEXT + +BEGIN_SOLUTION +\[ $dualtableau1_string \] +$PAR +Setting the parameters \(y_1, y_2\) to zero is not a feasible solution. +Options for finding a first feasible solution include guessing (not a bad choice +for a problem this small), using prior knowledge (e.g. a known optimal solution to +a similar problem), creating an auxiliary problem to find a feasible point, +and the "shortcut method" which is an accelerated version of the "auxiliary method". +END_SOLUTION + + +# ######################################################################################################### +# # Begin section 6 -- solve the tableau using simplex method +# # section 6 (13matrix, 14essay) +# # First we have to find a feasible solution -- phase 1 +# # We'll use the auxiliary method -- adding an extra slack variable. +# ######################################################################################################### + + +Context($context); +$dualtableau2aux = [[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1,0,0,0,-$bill_profit], + [-1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, 0,1,0,0,-$steve_profit], + [1,0,0,0,0,0,0,0,1,0]]; + +# [0,$money_total,$time_total,1,1,0,0,1,0,0]]; +$dualmatrix2aux = Matrix($dualtableau2aux); +$original_dual_aux_matrix=$dualmatrix2aux; +$dualtoplabelsphase1 = [qw(y0 y1 y2 y3 y4 y5 y6 v z b)]; +$dualmatrix2aux->{top_labels}= $dualtoplabelsphase1; +$dual_constraint1phase1 = Compute("-${bill_money_commitment}y1 -${bill_time_commitment}y2 -y3 <=-$bill_profit+y0"); +$dual_constraint2phase1 = Compute("-${steve_money_commitment}y1 -${steve_time_commitment}y2 -y4 <=-$steve_profit+y0"); +$dualconstraintsphase1 = List($dual_constraint1phase1,$dual_constraint2phase1); +$dualpivot2aux=Point(1,1); +$dualbasis2aux=Set(6,7); +######################################################################################################### +# section 6 answers (15essay 16matrix) +######################################################################################################### +Context()->texStrings; + +BEGIN_TEXT +Using the auxiliary method write the LOP for finding +the first feasible point. Write +the new constraints to be used by auxiliary method. $PAR +\{ ANS($dualconstraintsphase1->cmp),ans_box(3,80) \} +$PAR +We will continue to use the convention that we are maximizing +the auxiliary function so there will be some minus signs that need to be taken into account. +We'll let \(z = -y_0\) and try to maximize \(z\). +$PAR +Now construct the tableau you'll use for the auxiliary method. +Add a first column for the extra slack variable and a next to the last row +for the new objective function value \(z\) and write the new tableau. +The next to the last row will hold the objective function for \(z\). +The last row will hold the original +dual objective function which we'll just carry along. +$PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix2)\} +\{ANS($dualmatrix2aux->cmp()), $dualmatrix2aux->ans_array \} +$PAR +END_TEXT +######################################################################################################### +# debugging display code +######################################################################################################### +BEGIN_TEXT + +Display dualtableau2aux +\($dualtableau2aux\), matrix \($dualmatrix2aux\) +dualbasis \($dualbasis2aux\) and \($dualpivot2aux\) +$PAR +\{display_tableau_state($original_dual_aux_matrix,$dualtableau2aux, $dualmatrix2aux, $dualbasis2aux, $dualpivot2aux)\} +$PAR +END_TEXT + +Context()->normalStrings; + +Context()->texStrings; +BEGIN_SOLUTION +New objective: Maximize \(z = -y_0 \) subject to +\[\begin{align} +$dual_constraint1phase1&\\ +$dual_constraint2phase1 &\\ +\end{align} +\] with all variables non-negative. +$PAR +If \(y_0\) is large enough this always has a feasible solution with +\(y_1=y_2=0\). + +\[ \{lp_display_mm($dualtableau2aux) \} \] +$PAR +END_SOLUTION +Context()->normalStrings; + + + +# ######################################################################################################### +# # Section 7 final comparisons of primary and dual results. +# ######################################################################################################### + $dualpivot3aux=$dualpivot2aux=Point(1,1); #already set + ($dualtableau3aux,$dualbasis3aux, $dualstatevars3aux) = lp_basis_pivot($dualtableau2aux,$dualbasis2aux,$dualpivot2aux); + $B3 = matrix_from_submatrix($original_dual_aux_matrix,rows=>[1..$dual_row_size],columns=>($dualbasis3aux->data)); + $dualmatrix3aux = ($B3->det)*Matrix($dualtableau3aux); + $dualmatrix3aux->{top_labels}=$dualtoplabelsphase1; + $dualstate3aux=Matrix([[@$dualstatevars3aux]]); + + +# $dualpivot4aux = Point("(1,1)"); +# $dualtableau3aux = lp_pivot($dualtableau2,0,0); +# $dualmatrix3 = Matrix($dualtableau3aux); +# $dualmatrix3->{top_labels}=$dualtoplabelsphase1; +# $dualtableau3_string = lp_display_mm($dualtableau3aux); +# $z_initial=-$bill_profit; # a hint +# # this is the initial value of the z=-y0 variable we are maximizing? + +# now start to maximize z +$dualpivot4aux = Point("(2,2)"); +$dualtableau4aux = lp_pivot($dualtableau3aux,1,1); +$denom1 = 2000; +$dualmatrix4aux = $denom1*Matrix($dualtableau4aux); +$dualmatrix4aux->{top_labels}=$dualtoplabelsphase1; +$dualtableau4aux_string = lp_display_mm($dualmatrix4aux); + +$dualpivot5aux= Point("(1,3)"); +$dualtableau5aux = lp_pivot($dualtableau4aux,0,2); + +## FIXME entries close to zero in a matrix are not being compared properly. +## FIXME tolerance and tolType are being ignored + +#$dualtableau5->[2][3]=0; +#$dualtableau5->[2][5]=0; + +$denom2 = 1.3E6; +$dualmatrix5aux = $denom2*Matrix($dualtableau5aux); +$dualmatrix5aux->{top_labels}=$dualtoplabelsphase1; +$dualtableau5_string = lp_display_mm($dualmatrix5aux); +$state5 = Matrix([[ 0.423, 6.4615, 0, 0, -6415.38]]); +$popup5 = PopUp(["?","Yes", "No"], "No"); + +$dualpivot5aux= Point("(1,3)"); +$dualtableau5aux = lp_pivot($dualtableau4aux,0,2); + +## FIXME entries close to zero in a matrix are not being compared properly. +## FIXME tolerance and tolType are being ignored + +#$dualtableau5->[2][3]=0; +#$dualtableau5->[2][5]=0; + +$denom2 = 1.3E6; +$dualmatrix5aux = $denom2*Matrix($dualtableau5aux); +$dualmatrix5aux->{top_labels}=$dualtoplabelsphase1; +$dualtableau5_string = lp_display_mm($dualmatrix5aux); + +######################################################################################################### +# debugging display code +######################################################################################################### +BEGIN_TEXT +$PAR +(1,1) Debugging display code dualtableau3aux. This is the matrix +after pivoting on (1,1) \($dualpivot3aux\) to basis $dualbasis3aux . +newbasis: \{$newdualbasis3aux=Set(1,7)\} $newdualbasis3aux +$PAR + +\{display_tableau_state($original_dual_aux_matrix, + $dualtableau3aux, $dualmatrix3aux, + $newdualbasis3aux, $dualpivot3aux)\} +$PAR + +(2,2) Debugging display code dualtableau4aux. +This is the matrix after pivoting on $dualpivot4aux and basis $dualbasis3aux +and basis $dualbasis3aux . +newbasis: \{$newdualbasis4aux=Set(1,2)\} $newdualbasis4aux +$PAR +\{display_tableau_state($original_dual_aux_matrix, + $dualtableau4aux, $dualmatrix4aux, + $newdualbasis4aux, $dualpivot4aux)\} +$PAR + +(1,3) Debugging display code dualtableau5aux. +This is the matrix after pivoting on $dualpivot5aux and basis $dualbasis3aux . +newbasis: \{$newdualbasis5aux=Set(2,3)\} +$PAR +\{display_tableau_state($original_dual_aux_matrix, + $dualtableau5aux, $dualmatrix5aux, + $newdualbasis5aux, $dualpivot5aux)\} +$PAR +END_TEXT +# ######################################################################################################### +# # Section 7 answers 17,18matrix,19,20matrix,21,22matrix,23, 24,25essay +# ######################################################################################################### + + + +Context()->texStrings; + +BEGIN_TEXT +The first pivot, which will make the right hand side entries positive is +\{ANS($dualpivot3aux->cmp),$dualpivot3aux->ans_rule \}. $PAR +Recall that in the case of tie we take entry with the least index +(i.e. left most or upper most. ) +$PAR + + +The resulting tableau is $PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix3aux)\} +\{ANS($dualmatrix3aux->cmp), $dualmatrix3aux->ans_array() \} + +The value of \(z=-y_0\) is \($z_initial\). +$PAR +The next pivot, following the simplex method to maximize \(z\), is +\{ANS($dualpivot4aux->cmp), $dualpivot4aux->ans_rule\}. +$PAR +(Many of the following answers have lots of zeros. You can use the shortcut +1E3 to stand for \(1\times 10^3\). ) +$PAR +Notice that because of the zero on the right hand side none of the state variables change. We had three +hyperplanes intersecting at a point and we have changed our mind about which +of those three we consider basic. The new tableau is +$PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix4aux)\} +\{ANS($dualmatrix4aux->cmp), $dualmatrix4aux->ans_array() \} +$PAR +The next pivot \{ ANS($dualpivot5aux->cmp), $dualpivot5aux->ans_rule \} leads to $PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix5aux)\} +\{#FIXME -- these tolerance leveals are being ignored \} +\{ANS($dualmatrix5aux->cmp()), $dualmatrix5aux->ans_array() \} + +$PAR and we notice that now \(z=-y_0=0\) so we have found a basic feasible solution to our +original dual problem. The variables \(y_1,y_2,y_3,y_4,v\) for this +solution are +$PAR + +\{ANS($state5->cmp),$state5->ans_array\} +$PAR +Do we need to continue to optimize the value for \(v\)? +\{ANS($popup5->cmp), $popup5->menu()\} Why? $PAR +\{ANS(essay_cmp), essay_box(3,80)\} +$PAR +Compare this answer \(v^*= -w^*\) to the dual problem to the optimal value \(P^*\)for the primary problem. +The problems' goals were to maximize \(P\) and to minimize \(w\). +$PAR + +END_TEXT + +BEGIN_SOLUTION +The first pivot is \($dualpivot3aux\) which makes the right hand side entries positive. +$PAR +\[$dualtableau3_string\] +$PAR +The next pivot is \($dualpivot4aux\). It follows the simplex rule of choosing the row +with the most restrictive ratio -- in this case zero. The result is simply to +choose a new representation of the same point -- no change in state takes place. +$PAR\[$dualtableau4_string\] +$PAR +The final pivot (as it turns out) is \($dualpivot5aux\). At this point we are done with +phase 1 because the value of \(z=-y_0=0\) so we have maximized \(z\) and +minimized \(y_0\) to \(0\) +$PAR\[$dualtableau5_string\] +$PAR +At this point we notice that the coefficients in the last row are such that +we cannot increase the value of \(v\) any further. \(v\) is at its maximum +and \(w=-v\) is at its minimum. The minimum value of \(w=64?\) which is the +same as the maximum profit that we found in the first example. +END_SOLUTION +TEXT($END_ONE_COLUMN); +ENDDOCUMENT(); + diff --git a/t/matrix_tableau_tests/tableau_test5.pg b/t/matrix_tableau_tests/tableau_test5.pg new file mode 100644 index 0000000000..51648e8c7a --- /dev/null +++ b/t/matrix_tableau_tests/tableau_test5.pg @@ -0,0 +1,505 @@ +DOCUMENT(); +loadMacros( +"PGstandard.pl", +"MathObjects.pl", +"parserPopUp.pl", +"unionLists.pl", +"MatrixReduce.pl", +#"AppletObjects.pl", +"PGessaymacros.pl", +"PGmatrixmacros.pl", +"LinearProgramming.pl", +"parserLinearInequality.pl", +"quickMatrixEntry.pl", +#"scaffold.pl", +"tableau.pl", +#"gage_matrix_ops.pl", +"PGinfo.pl", +"source.pl", +"PGcourse.pl", +); + +TEXT(beginproblem()); +TEXT($BEGIN_ONE_COLUMN); +$showPartialCorrectAnswers = 1; + +INITIALIZE_QUICK_MATRIX_ENTRY(); + + +############################################################## +# problem data +############################################################## + +# this can be changed (slightly at least) without affecting the behavior of the problem +# The choice of pivots is not yet automatic However + +# Your resources: +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + +#Hack to prevent domain conflict in answer. +Context()->variables->add(p1=>'Real',p2=>'Real'); +$objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); +# why can't the formula be defined within context "linearInequality"? + +Context("LinearInequality"); +Context()->variables->add(p1=>'Real',p2=>'Real'); +Context()->strings->add("Essay Answer" =>{}); +Context()->strings->add('Minimize'=>{},'Maximize'=>{}, "?"=>{}); +Context()->strings->add('Yes'=>{},'No'=>{}); +Context()->flags->set({tolType=>"absolute",tolerance=>.001}); +Context()->flags->set( + zeroLevel=>0.001, + zeroLevelTol=>.001 +); +# $objfun1 = Formula("${bill_profit}p1 + ${steve_profit}p2"); + +our $context=Context(); + +############################################################## +# +# Setup +# +# + + + +############################################################################## +# utility subroutine for checking your answers +############################################################################## + +sub update_tableau { + my $original_matrix = shift; + my $basis = shift; + # $basis is a Set object + # last row is numbered $basis_size+1; + my $basis_size=@{$basis->data}; + my $basis_matrix = $original_matrix->submatrix(rows=>[1..$basis_size],columns=>$basis->data); + + my $obj_coeff = - ($original_matrix->row($basis_size+1)); + #obj function coeff corresponding to basis variables/columns + $cB = Matrix($obj_coeff->column_slice($basis->data)); + # get A | S |b portion of tableau and multiply by det(B)B^-1 + my $current_matrix= ($basis_matrix->det)* + ($basis_matrix->inverse)* + $original_matrix->row_slice(1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + return $current_tableau; +} +# the full tableau has one more row -- the objective function +sub display_tableau_state { + my ($original_matrix,$tableau, $matrix, $basis,$pivot)= @_; + $basis = $basis->sort; + my $basis_size = @{$basis->data} ; + my $basis_matrix = matrix_from_submatrix($original_matrix, rows=>[1..$basis_size],columns=>$basis->data); + my $normalized_tableau = $matrix; + my $reduced_matrix = matrix_from_matrix_rows($matrix,1..$basis_size); + warn "basis matrix is $basis_matrix, $reduced_matrix"; + my $normalized_reduced_matrix = ($basis_matrix->det)*$reduced_matrix; + # calculate current tableau by multiplying by the inverse of the basis + my $obj_coeff = - $original_matrix->row($basis_size+1); + my $cB = Matrix($obj_coeff->column_slice($basis->data)); + my $current_matrix = ($basis_matrix->det)*($basis_matrix->inverse)*matrix_from_matrix_rows($original_matrix,1..$basis_size); + my $new_obj_row = ($cB*$current_matrix)->row(1) - ($basis_matrix->det)*($obj_coeff); + my $current_tableau = Matrix(@{$current_matrix->extract_rows},Matrix($new_obj_row)); + $current_tableau = update_tableau($original_matrix, $basis); + + #DEBUG_MESSAGE("get state variables for $matrix using basis $basis"); + @statevars1 = get_tableau_variable_values($matrix, $basis); + # get z value + $statevars1 = ~~@statevars1; + $state = Matrix([[@statevars1]]); + + + return " + + pivot: \($pivot\) basis: \( $basis\) $PAR basis matrix \(". + display_matrix_mm($basis_matrix)."\) $PAR + tableau and normalized current tableau $PAR + \(" . lp_display_mm($tableau) . "\)\(".lp_display_mm($normalized_tableau). " \) $PAR + original matrix and current matrix calculated by Binverse*original $PAR + \(". + display_matrix_mm(matrix_from_matrix_rows($original_matrix,1..$basis_size)). + "\) \(". + display_matrix_mm($current_matrix). + "\) $PAR + objective coefficients = \( $obj_coeff \) $PAR + basis objective_coefficients = \($cB\) $PAR + new tableau \[". + lp_display_mm($current_tableau). + "\] $PAR + state: \($state)\) $PAR + " +} + + +# ######################################################################################################### +# # Section 4 dual problem +# ######################################################################################################### + + +Context()->variables->add(y1=>'Real',y2=>'Real', y3=>'Real', y4=>'Real', w=>'Real',y0=>'Real'); +$dual_constraint1 = Compute("${bill_money_commitment}y1 +${bill_time_commitment}y2 +y3 >=$bill_profit"); +$dual_constraint2 = Compute("${steve_money_commitment}y1 +${steve_time_commitment}y2 +y4 >=$steve_profit"); +$dual_constraints= List($dual_constraint1, $dual_constraint2); + +$dual_objfun = Formula("${money_total}y1 +${time_total}y2 + y3 +y4"); +$popupmaxmin = PopUp(["?","Maximize", "Minimize"], "Minimize"); + +Context()->texStrings; + +BEGIN_TEXT +Construct the dual problem for the linear optimization problem above. The first goal is to calculate +an upper bound for the possible profit in the LOP using linear combinations of the inequality constraints. +Then formulate the search for the best of these upper bounds in such a way that it +becomes a new LOP -- the dual problem. Use variables \(y_1,y_2,y_3,y_4\) to create linear combinations +of the constraints on money, time and the total probabilities \(p_1\) and \(p_2\) respectively, +and create a linear function \(w = Ay_1 + By_2 +Cy_3 +Dy_4\) which guarantees that \(w\) will +be larger than any profit one could make. What constraints must \(y_1\dots y_4\) satisfy? +$PAR +\{ANS($dual_constraints->cmp->withPostFilter( + linebreak_at_commas() +)),ans_box(4,80)\} +$PAR The objective function would be: + \(w = \) \{ANS($dual_objfun->cmp), $dual_objfun->ans_rule(50)\} +$PAR +To get the best possible estimate for the profit in the original problem +you would want to \{ANS($popupmaxmin->cmp),$popupmaxmin->menu\} \(w\) +Why? +$PAR +\{ANS(essay_cmp()), essay_box(5,80)\} +END_TEXT + +BEGIN_SOLUTION +The profit is less than the minimum of \( w = $dual_objfun \) subject to + +\[\begin{aligned} +$dual_constraint1 & \\ +$dual_constraint2 & +\end{aligned} +\] +The \(y_i\) are the coefficients used to add up the constraint inequalities of the primary problem +so as to estimate an upper bound for the profit. If the \(y_i\) are positive values satisfying +the constraints then \(w\) will be greater than or equal to the profit. We get the most precise +estimate by finding values of \(y_i\), satisfying the constraints, which give the smallest +possible value for \( w\). + +END_SOLUTION +Context()->normalStrings; + + +# ######################################################################################################### +# # Section 5 set up dual constraints +# ######################################################################################################### + +Context($context); +$dualtableau1 = [[-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1, 0, 0,-$bill_profit], + [-$steve_money_commitment, -$steve_time_commitment,0,-1, 0, 1, 0,-$steve_profit], + [${money_total},${time_total},1, 1, 0,0, 1,0]]; +$dualmatrix1 = Matrix($dualtableau1); +$dualtoplabels = [qw(y1 y2 y3 y4 y5 y6 v b)]; +$dualmatrix1->{top_labels}=$dualtoplabels; +$dualtableau1_string = lp_display_mm($dualtableau1,top_labels=>$dualtoplabels); + +$popup = PopUp(["?","Yes", "No"], "No"); + +$dual_row_size =2; # figure these out automatically +$dual_col_size =$col_size; +BEGIN_TEXT +dualmatrix1 size \{join(",", $dualmatrix1->dimensions)\} $BR +$dual_row_size, $dual_col_size +END_TEXT +Context()->texStrings; + +BEGIN_TEXT +$PAR +For this first effort let's rewrite the equations so that they are in "standard" form, meaning +that the inequalities are all less than or equal to and that the goal is to maximize \(v = -w\). +This will probably mean that you will have to rewrite the natural way in which you set up +the problem above. Some coefficients will change sign. +We will still be using the convention that +the coefficient of \(v\) will be \(1\). + +$PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix1)\} +\{ANS($dualmatrix1->cmp()), $dualmatrix1->ans_array \} +$PAR +(This tableau is for maximizing \(v\), aka \(-w\) -- there are two minus sign switches. ) +$PAR +Is there a natural feasible solution to this problem? In other words does setting the +problem parameters equal to 0 provide a feasible solution? +\{ANS($popup->cmp),$popup->menu \} +$PAR +Explain why or why not and what are your options for getting started. +\{ ANS(essay_cmp()), essay_box(10,80) \} +$PAR +END_TEXT + +BEGIN_SOLUTION +\[ $dualtableau1_string \] +$PAR +Setting the parameters \(y_1, y_2\) to zero is not a feasible solution. +Options for finding a first feasible solution include guessing (not a bad choice +for a problem this small), using prior knowledge (e.g. a known optimal solution to +a similar problem), creating an auxiliary problem to find a feasible point, +and the "shortcut method" which is an accelerated version of the "auxiliary method". +END_SOLUTION + + +# ######################################################################################################### +# # Begin section 6 -- solve the tableau using simplex method +# # section 6 (13matrix, 14essay) +# # First we have to find a feasible solution -- phase 1 +# # We'll use the auxiliary method -- adding an extra slack variable. +# ######################################################################################################### + + +Context($context); +$dualtableau2aux = [[-1,-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1,0,0,0,-$bill_profit], + [-1,-$steve_money_commitment,-$steve_time_commitment, 0, -1, 0,1,0,0,-$steve_profit], + [1,0,0,0,0,0,0,0,1,0]]; + +# [0,$money_total,$time_total,1,1,0,0,1,0,0]]; +$dualmatrix2aux = Matrix($dualtableau2aux); +$original_dual_aux_matrix=$dualmatrix2aux; +$dualtoplabelsphase1 = [qw(y0 y1 y2 y3 y4 y5 y6 v z b)]; +$dualmatrix2aux->{top_labels}= $dualtoplabelsphase1; +$dual_constraint1phase1 = Compute("-${bill_money_commitment}y1 -${bill_time_commitment}y2 -y3 <=-$bill_profit+y0"); +$dual_constraint2phase1 = Compute("-${steve_money_commitment}y1 -${steve_time_commitment}y2 -y4 <=-$steve_profit+y0"); +$dualconstraintsphase1 = List($dual_constraint1phase1,$dual_constraint2phase1); +$dualpivot2aux=Point(1,1); +$dualbasis2aux=Set(6,7); +######################################################################################################### +# section 6 answers (15essay 16matrix) +######################################################################################################### +Context()->texStrings; + +BEGIN_TEXT +Using the auxiliary method write the LOP for finding +the first feasible point. Write +the new constraints to be used by auxiliary method. $PAR +\{ ANS($dualconstraintsphase1->cmp),ans_box(3,80) \} +$PAR +We will continue to use the convention that we are maximizing +the auxiliary function so there will be some minus signs that need to be taken into account. +We'll let \(z = -y_0\) and try to maximize \(z\). +$PAR +Now construct the tableau you'll use for the auxiliary method. +Add a first column for the extra slack variable and a next to the last row +for the new objective function value \(z\) and write the new tableau. +The next to the last row will hold the objective function for \(z\). +The last row will hold the original +dual objective function which we'll just carry along. +$PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix2)\} +\{ANS($dualmatrix2aux->cmp()), $dualmatrix2aux->ans_array \} +$PAR +END_TEXT +######################################################################################################### +# debugging display code +######################################################################################################### +BEGIN_TEXT + +Display dualtableau2aux +\($dualtableau2aux\), matrix \($dualmatrix2aux\) +dualbasis \($dualbasis2aux\) and \($dualpivot2aux\) +$PAR +\{display_tableau_state($original_dual_aux_matrix,$dualtableau2aux, $dualmatrix2aux, $dualbasis2aux, $dualpivot2aux)\} +$PAR +END_TEXT + +Context()->normalStrings; + +Context()->texStrings; +BEGIN_SOLUTION +New objective: Maximize \(z = -y_0 \) subject to +\[\begin{align} +$dual_constraint1phase1&\\ +$dual_constraint2phase1 &\\ +\end{align} +\] with all variables non-negative. +$PAR +If \(y_0\) is large enough this always has a feasible solution with +\(y_1=y_2=0\). + +\[ \{lp_display_mm($dualtableau2aux) \} \] +$PAR +END_SOLUTION +Context()->normalStrings; + + + +# ######################################################################################################### +# # Section 7 final comparisons of primary and dual results. +# ######################################################################################################### + $dualpivot3aux=$dualpivot2aux=Point(1,1); #already set + ($dualtableau3aux,$dualbasis3aux, $dualstatevars3aux) = lp_basis_pivot($dualtableau2aux,$dualbasis2aux,$dualpivot2aux); + $B3 = matrix_from_submatrix($original_dual_aux_matrix,rows=>[1..$dual_row_size],columns=>($dualbasis3aux->data)); + $dualmatrix3aux = ($B3->det)*Matrix($dualtableau3aux); + $dualmatrix3aux->{top_labels}=$dualtoplabelsphase1; + $dualstate3aux=Matrix([[@$dualstatevars3aux]]); + + +# $dualpivot4aux = Point("(1,1)"); +# $dualtableau3aux = lp_pivot($dualtableau2,0,0); +# $dualmatrix3 = Matrix($dualtableau3aux); +# $dualmatrix3->{top_labels}=$dualtoplabelsphase1; +# $dualtableau3_string = lp_display_mm($dualtableau3aux); +# $z_initial=-$bill_profit; # a hint +# # this is the initial value of the z=-y0 variable we are maximizing? + +# now start to maximize z +$dualpivot4aux = Point("(2,2)"); +$dualtableau4aux = lp_pivot($dualtableau3aux,1,1); +$denom1 = 2000; +$dualmatrix4aux = $denom1*Matrix($dualtableau4aux); +$dualmatrix4aux->{top_labels}=$dualtoplabelsphase1; +$dualtableau4aux_string = lp_display_mm($dualmatrix4aux); + +$dualpivot5aux= Point("(1,3)"); +$dualtableau5aux = lp_pivot($dualtableau4aux,0,2); +$denom2 = 1.3E6; +$dualmatrix5aux = $denom2*Matrix($dualtableau5aux); +$dualmatrix5aux->{top_labels}=$dualtoplabelsphase1; +$dualtableau5_string = lp_display_mm($dualmatrix5aux); +$state5 = Matrix([[ 0.423, 6.4615, 0, 0, -6415.38]]); +$popup5 = PopUp(["?","Yes", "No"], "No"); + +$dualpivot5aux= Point("(1,3)"); +$dualtableau5aux = lp_pivot($dualtableau4aux,0,2); + +## FIXME entries close to zero in a matrix are not being compared properly. +## FIXME tolerance and tolType are being ignored + +#$dualtableau5->[2][3]=0; +#$dualtableau5->[2][5]=0; + +$denom2 = 1.3E6; +$dualmatrix5aux = $denom2*Matrix($dualtableau5aux); +$dualmatrix5aux->{top_labels}=$dualtoplabelsphase1; +$dualtableau5_string = lp_display_mm($dualmatrix5aux); + +######################################################################################################### +# debugging display code +######################################################################################################### +BEGIN_TEXT +$PAR +(1,1) Debugging display code dualtableau3aux. This is the matrix +after pivoting on (1,1) \($dualpivot3aux\) to basis $dualbasis3aux . +newbasis: \{$newdualbasis3aux=Set(1,7)\} $newdualbasis3aux +$PAR + +\{display_tableau_state($original_dual_aux_matrix, + $dualtableau3aux, $dualmatrix3aux, + $newdualbasis3aux, $dualpivot3aux)\} +$PAR + +(2,2) Debugging display code dualtableau4aux. +This is the matrix after pivoting on $dualpivot4aux and basis $dualbasis3aux +and basis $dualbasis3aux . +newbasis: \{$newdualbasis4aux=Set(1,2)\} $newdualbasis4aux +$PAR +\{display_tableau_state($original_dual_aux_matrix, + $dualtableau4aux, $dualmatrix4aux, + $newdualbasis4aux, $dualpivot4aux)\} +$PAR + +(1,3) Debugging display code dualtableau5aux. +This is the matrix after pivoting on $dualpivot5aux and basis $dualbasis3aux . +newbasis: \{$newdualbasis5aux=Set(2,3)\} +$PAR +\{display_tableau_state($original_dual_aux_matrix, + $dualtableau5aux, $dualmatrix5aux, + $newdualbasis5aux, $dualpivot5aux)\} +$PAR +END_TEXT +# ######################################################################################################### +# # Section 7 answers 17,18matrix,19,20matrix,21,22matrix,23, 24,25essay +# ######################################################################################################### + + + +Context()->texStrings; + +BEGIN_TEXT +The first pivot, which will make the right hand side entries positive is +\{ANS($dualpivot3aux->cmp),$dualpivot3aux->ans_rule \}. $PAR +Recall that in the case of tie we take entry with the least index +(i.e. left most or upper most. ) +$PAR + + +The resulting tableau is $PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix3aux)\} +\{ANS($dualmatrix3aux->cmp), $dualmatrix3aux->ans_array() \} + +The value of \(z=-y_0\) is \($z_initial\). +$PAR +The next pivot, following the simplex method to maximize \(z\), is +\{ANS($dualpivot4aux->cmp), $dualpivot4aux->ans_rule\}. +$PAR +(Many of the following answers have lots of zeros. You can use the shortcut +1E3 to stand for \(1\times 10^3\). ) +$PAR +Notice that because of the zero on the right hand side none of the state variables change. We had three +hyperplanes intersecting at a point and we have changed our mind about which +of those three we consider basic. The new tableau is +$PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix4aux)\} +\{ANS($dualmatrix4aux->cmp), $dualmatrix4aux->ans_array() \} +$PAR +The next pivot \{ ANS($dualpivot5aux->cmp), $dualpivot5aux->ans_rule \} leads to $PAR +\{MATRIX_ENTRY_BUTTON($dualmatrix5aux)\} +\{#FIXME -- these tolerance leveals are being ignored \} +\{ANS($dualmatrix5aux->cmp()), $dualmatrix5aux->ans_array() \} + +$PAR and we notice that now \(z=-y_0=0\) so we have found a basic feasible solution to our +original dual problem. The variables \(y_1,y_2,y_3,y_4,v\) for this +solution are +$PAR + +\{ANS($state5->cmp),$state5->ans_array\} +$PAR +Do we need to continue to optimize the value for \(v\)? +\{ANS($popup5->cmp), $popup5->menu()\} Why? $PAR +\{ANS(essay_cmp), essay_box(3,80)\} +$PAR +Compare this answer \(v^*= -w^*\) to the dual problem to the optimal value \(P^*\)for the primary problem. +The problems' goals were to maximize \(P\) and to minimize \(w\). +$PAR + +END_TEXT + +BEGIN_SOLUTION +The first pivot is \($dualpivot3aux\) which makes the right hand side entries positive. +$PAR +\[$dualtableau3_string\] +$PAR +The next pivot is \($dualpivot4aux\). It follows the simplex rule of choosing the row +with the most restrictive ratio -- in this case zero. The result is simply to +choose a new representation of the same point -- no change in state takes place. +$PAR\[$dualtableau4_string\] +$PAR +The final pivot (as it turns out) is \($dualpivot5aux\). At this point we are done with +phase 1 because the value of \(z=-y_0=0\) so we have maximized \(z\) and +minimized \(y_0\) to \(0\) +$PAR\[$dualtableau5_string\] +$PAR +At this point we notice that the coefficients in the last row are such that +we cannot increase the value of \(v\) any further. \(v\) is at its maximum +and \(w=-v\) is at its minimum. The minimum value of \(w=64?\) which is the +same as the maximum profit that we found in the first example. +END_SOLUTION +TEXT($END_ONE_COLUMN); +ENDDOCUMENT(); + diff --git a/t/matrix_tableau_tests/test_scalar_mult.pg b/t/matrix_tableau_tests/test_scalar_mult.pg new file mode 100644 index 0000000000..0be93bcc4b --- /dev/null +++ b/t/matrix_tableau_tests/test_scalar_mult.pg @@ -0,0 +1,71 @@ +##DESCRIPTION +# +# +# tests multiplying a matrix by a large scalar inside BEGIN_TEXT/END_TEXT +# if the scalar is large enough this fails in the presence of texStrings because +# the tex description of $c*$m becomes 1.3\times 10^{6}*\left[\begin{array}{cc} +# 1 &2\cr +# 3 &4 +# \end{array}\right] +# which can't be +##ENDDESCRIPTION + + +DOCUMENT(); # This should be the first executable line in the problem. + +loadMacros( + "PGstandard.pl", # Standard macros for PG language + "MathObjects.pl", + #"source.pl", # used to display problem source button + "PGcourse.pl", # Customization file for the course +); + +TEXT(beginproblem()); +$showPartialCorrectAnswers = 1; + +############################################################## +# +# Setup +# +# +Context("Matrix"); + +############################################################## +# +# Text +# +# + + +# in the current develop branch (10/14/2017) there should be no errors generated. + + +$m = Matrix([[0,1,0,0,0,0,0,0,0,0.423077],[0.00153846,0,1,-0.00230769,0.00384615,0.00230769,-0.00384615,0,0,6.46154]]); +$c = Real(1300000); # this gives "can't convert a Real Number to a Matrix" when executed inside texStrings mode; +# $c = Real(130000); # this works +#$c = Real(1300000)->value; # this also works +#$m = Matrix([[1,2],[3,4]]); + +$string = $c*$m; # this does not cause errors. + +Context()->texStrings; +$string = $c*$m; # this causes errors. +BEGIN_TEXT + +\[ $string \] $PAR + +\[\{ $c*$m \} \] + +END_TEXT +Context()->normalStrings; + + +############################################################## +# +# Answers +# +# + + + +ENDDOCUMENT(); # This should be the last executable line in the problem. \ No newline at end of file diff --git a/t/matrix_tableau_tests/test_tableau_more.pl b/t/matrix_tableau_tests/test_tableau_more.pl new file mode 100755 index 0000000000..3a1873a312 --- /dev/null +++ b/t/matrix_tableau_tests/test_tableau_more.pl @@ -0,0 +1,69 @@ +#!/usr/bin/perl -w + +use lib "/Volumes/WW_test/opt/webwork/pg_2014/lib"; +use lib "/Volumes/WW_test/opt/webwork/pg_2014/macros"; +use lib "/Volumes/WW_test/opt/webwork/webwork2/lib"; + +use Test::More; +use Parser; +use Value; +require "tableau.pl"; +require "Value.pl"; + + + +$money_total = 6000; +$time_total = 600; + +# Bill +$bill_money_commitment = 5000; #dollars +$bill_time_commitment = 400; # hours +$bill_profit = 4700; +# Steve +$steve_money_commitment = 3000; +$steve_time_commitment = 500; +$steve_profit = 4500; + + + +#### problem starts here: + +# need error checking to make sure that tableau->new checks +# that inputs are matrices +$ra_matrix = [[-$bill_money_commitment,-$bill_time_commitment, -1, 0, 1,0,0,-$bill_profit], + [-$steve_money_commitment,-$steve_time_commitment, 0, -1, 0,1,0,-$steve_profit], + [-$money_total,-$time_total,0,0, 0,0, 1,0]]; + +$A = Value::Matrix->new([[-$bill_money_commitment,-$bill_time_commitment, -1, 0], + [ -$steve_money_commitment,-$steve_time_commitment, 0, -1 ]]); +$b = Value::Vector->new([-$bill_profit,-$steve_profit]); # need vertical vector +$c = Value::Vector->new([$money_total,$time_total,0,0]); + +$tableau1 = Tableau->new(A=>$A, b=>$b, c=>$c); + +ok (1==1, "trivial first test"); +ok (defined($tableau1), 'tableau has been defined and loaded'); +is ($tableau1->{m}, 2, 'number of constraints is 2'); +is ($tableau1->{n}, 4, 'number of variables is 4'); +is_deeply ( [$tableau1->{m},$tableau1->{n}], [$tableau1->{A}->dimensions], '{m},{n} match dimensions of A'); +is_deeply ($tableau1->{A}, $A, 'constraint matrix'); +is_deeply ($tableau1->{b}, Matrix([$b])->transpose, 'constraint constants is m by 1 matrix'); +is_deeply ($tableau1->{c}, $c, 'objective function constants'); +my $test_constraint_matrix = Matrix($ra_matrix->[0],$ra_matrix->[1]); +is_deeply ($tableau1->{current_constraint_matrix}, $test_constraint_matrix, + 'initialization of current_constraint_matrix'); +is_deeply ($tableau1->{current_b}, $tableau1->{b}, + 'initialization of current_b'); +my $obj_row_test = [ ((-$c)->value, 0,0,1,0) ]; +is_deeply ($tableau1->objective_row, $obj_row_test, + 'initialization of $tableau->{obj_row}'); +is_deeply( ref($tableau1->objective_row), 'ARRAY', '->objective_row has type ARRAY'); +is_deeply( ref($tableau1->{obj_row}), 'Value::Matrix', '->{obj_row} has type Value::Matrix'); + +is_deeply( $tableau1->objective_row, [$tableau1->{obj_row}->value], 'access to {obj_row}'); +is_deeply($tableau1->current_tableau, Matrix($ra_matrix), 'entire tableau including obj coeff row'); +is(ref($tableau1->current_tableau), 'Value::Matrix', '-> current_tableau is Value::Matrix'); + +# check accessors? Use antlers? + +done_testing(); \ No newline at end of file