diff --git a/.gitignore b/.gitignore index b47981c..3ae937e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,39 +1,40 @@ -*.annot -*.cmo -*.cma -*.cmi -*.a -*.o -*.cmx -*.cmxs -*.cmxa -*.aux -*.bbl -*.blg -*.dvi -*.log -*.out -*.pdf -*.ps -*.gz -*~ -\#*\# -/.emacs.desktop -/.emacs.desktop.lock -.elc -auto-save-list -tramp -.\#* - -# Org-mode -.org-id-locations -*_archive -*.swp - -# packages -*.tar - -# program -simplex -matrixtop -simplextop +*.annot +*.cmo +*.cma +*.cmi +*.a +*.o +*.cmx +*.cmxs +*.cmxa +*.aux +*.bbl +*.blg +*.dvi +*.log +*.out +*.pdf +*.ps +*.gz +*~ +\#*\# +/.emacs.desktop +/.emacs.desktop.lock +.elc +auto-save-list +tramp +.\#* + +# Org-mode +.org-id-locations +*_archive +*.swp + +# packages +*.tar + +# program +maketop +simplex +matrixtop +simplextop diff --git a/DraftSpec_Annotated.tex b/DraftSpec_Annotated.tex index eafc62f..7f192d2 100644 --- a/DraftSpec_Annotated.tex +++ b/DraftSpec_Annotated.tex @@ -1,392 +1,392 @@ -\documentclass[letterpaper,11pt]{article} -% Change the header if you're going to change this. - -% Possible packages - uncomment to use -%\usepackage{amsmath} % Needed for math stuff -%\usepackage{lastpage} % Finds last page -%\usepackage{amssymb} % Needed for some math symbols -%\usepackage{graphicx} % Needed for graphics -\usepackage[usenames,dvipsnames]{xcolor} % Needed for graphics and color -%\usepackage{setspace} % Needed to set single/double spacing -\usepackage{hyperref} - -%Suggested by TeXworks -\usepackage[utf8]{inputenc} % set input encoding (not needed with XeLaTeX) -\usepackage{verbatim} % adds environment for commenting out blocks of text & for better verbatim - -% Sets the page margins to 1 inch each side -\usepackage[margin=1in]{geometry} - -% Uncomment this section if you wish to have a header. -\usepackage{fancyhdr} -\pagestyle{fancy} -\renewcommand{\headrulewidth}{0.5pt} % customise the layout... -\lhead{CS51 Draft Spec} \chead{} \rhead{May 5th, 2013} -\lfoot{} \cfoot{\thepage} \rfoot{} - -\frenchspacing -\setlength\parindent{0pt} -\setlength\parskip{2ex} - -\newcommand{\annot}[1]{\textbf{\color{BrickRed} [#1]}} - -\begin{document} -\title{CS 51 Final Project Functional Specification{\annot{Matrix Module and Simplex Algorithm}} ----{\annot{Annotated!}}} -\author{ -Luis Perez $|$ luisperez@college.harvard.edu \\ -Andy Shi $|$ andyshi@college.harvard.edu \\ -Zihao Wang $|$ zihaowang01@college.harvard.edu \\ -Ding Zhou $|$ dzhou@college.harvard.edu -} -\date{May 5, 2013} -\maketitle - -%Put your stuff here -{\annot{Note for the reader: This is the annotated version of our CS 51 -Final Project Functional Specifications. All annotations will appear in this color.}}. - -\section{Brief Overview} - -Our project seeks to implement an efficient matrix library for Ocaml. The end -goal is to utilize this matrix library in conjunction with the simplex algorithm -to find solution {\annot{Linear Programs}}. In order to accomplish this, -there will be multiple small checkpoints. We need to implement and efficient -matrix representation (so as to not use the naive implementation of lists of -lists){\annot{(Successfully accomplished. Our implementation uses array for efficienty)}}. - We also hope to provide a useful interface that allows the user to input -arbitrarily large matrices effectively (possible by mapping over the entries of -the matrix){\annot{The user can input matrices with a simple text file composed of comma -seperated values}}. One big challenge for the matrix library will be the implementation -of efficient algorithms for multiplication and row reduction. We expect the -multiplication to initially consists of the standard algorithm, but will -eventually grow into either block multiplication or Strassen {\annot{We were not -able to implement more sophisticated multiplication algorithms, but it seems that the current -algorithm works fairly well}}. We might even implement some way to chose between the -algorithms depending on the inputted matrices, but that would not an essential feature -{\annot{This feature can still easily be implemented, but there simply was no time. We -already knew this, though, by the time we wrote our Technical Specs}}. - -The row reduction algorithm will certainly be one of the most fundamental -algorithms in our matrix library implementation. From our current understanding -of simplex, it seems like it will heavily rely on row reduction {\annot{The -simplex algorithm ended up relying on a different type of row reduction, so we -couldn't just plug in a matrix into our row reduction algorithm. We had to do a -lot of extra work to get simplex. Though, many of the helper functions for row -reduction ended up being extremely useful for Simplex, so this was not as -pointless as some might think}}. The only algorithm we can come up with so far -is Gauss-Jordan, though we have found some algorithms which are more suited for -computation by computers with imprecise storage of floats {\annot{In fact, we -just had to be careful which element to choose as a pivot. Futhermore, we ended -up implementing arbitrary precision, via Ocaml's num library, so this became a -non-issue}}. - -The simplex part of the problem is the main goal, but if we can get down the -matrix module, then it appears as if the simplex algorithm should be -straightforward {\annot{This was definitely not the case}}. Eventually, we wish -to implement a simplex module that allows the input of arbitrary linear -equations with linear constraints, and then we parse this information into the -corresponding matrix and finally perform the simplex algorithm on that matrix in -order to return a solution to the specified problem.{\annot{This is fully functionalt. -Input can be made in the form of a text file with a linear program. The algorithm -runs and not only returns the right values for simple cases, but also checks for -cases that might initially appear to be unsolvable but aren't, and for cases which -are unbounded.}} - -Other tentative goals are to implement matrices that can accept either -arbitrarily large numbers (bignums) or arbitrarily precise floats -(\verb@int list * int stream@). It would also be interesting, though not -necessary, to implement multiple algorithms for each operation and choose from -them the one that is the best for a specific input data {\annot{We were able to -implement arbitrary precisions and arbitrarily large sizes with the OCaml Nums -library}}. - -\section{Feature List} -\begin{enumerate} - -\item Find a way to represent a matrix (arrays or list list or other -representation)---will likely use an array representation using the Ocaml array -module (\url{http://caml.inria.fr/pub/docs/manual-ocaml/libref/Array.html}) -{\annot{Implemented as specified}} - -\item Make a matrix from a list of lists (this is how the user will input the -data of the matrix). {\annot{Implemented as specified}} - -\item Scalar multiplication---as described here -(\url{http://www.purplemath.com/modules/mtrxmult.htm}) -{\annot{Implemented as specified}} - -\item Matrix addition---as described here -(\url{http://www.purplemath.com/modules/mtrxmult.htm}) -{\annot{Implemented as specified}} - -\item Standard Matrix multiplication -(\url{http://www.purplemath.com/modules/mtrxmult.htm}) -{\annot{Implemented as specified}} - -\item Row Reduction (Gauss-Jordan). We will use the following resource to help -us calculate row reduction when the entries of the matrix are stored with -limited precision -(\url{http://thejuniverse.org/PUBLIC/LinearAlgebra/LOLA/rowRed/var.html}) -{\annot{Implemented as specified}} - -\item Matrix inverse / Row reduction with LU Decomposition - Description can be -found here (\url{http://www.math.ust.hk/~macheng/math111/LU_Decomposition.pdf}) -{\annot{To find the inverse, we augmented the original matrix with an identity -matrix, then row-reduced the augmented matrix}} - -\item{\annot{Make a matrix from a text file. The user is now able to read in a -text file and create a matrix.}} - -\item{\annot{For Information of further functionality, take a look at MatrixI.ml.}} - -\item Simplex Algorithm (Described in Algorithms book, in addition to the -resources provided below) -{\annot{Implemented mostly as specified, with some modifications from the -Wikipedia article on simplex}} - -\item Matrix norm -{\annot{Did not implement since it was not necessary for Simplex.}} - -\item Matrix transpose -{\annot{Implemented as specified.}} - -\item Trace---This is a typical trace (sum of the entries on the diagonal) -{\annot{Implemented as specified.}} - -\item Determinant---found using algorithm described in Section 4.8 of Hubbard's -\emph{Vector Calculus, Linear Algebra, and Differential Forms: A Unified -Approach} (Math 23a/b textbook) -{\annot{Implemented by decomposing the matrix into upper and lower triangular -matrices, then multiplying the entries along the diagonals for both and -multiplying the two results}} - -\item Eigenvectors / Eigenvalues, calculated also by a method in Hubbard's book -(this might not be necessary to carry out the simplex algorithm and we might -abandon it if we can't figure out how to implement it in a timely manner) -{\annot{Did not implement since it was not necssary for Simplex.}} - -\end{enumerate} - -\section{Draft Technical Specification} - -\subsection{Matrices} - -Our matrix module will actually be a functor that takes in an -\verb@ORDERED_AND_OPERATIONAL@ module (from ps4 and moogle) so we can have -matrices of ints, floats, or floats as streams. -{\annot{This actually ended up being very helpful since at the end we we able -to use OCaml's Num Library in order to achieve arbitrary precision and size.}} - -The following type definition (a la ps4) will give us a way to talk about order -for different data types. - -\begin{verbatim} -type order = Equal | Less | Greater -\end{verbatim} - -The following signature will allow us to provide matrix functionality on a -variety of data types. For example, we could make an -\verb@ORDERED_AND_OPERATIONAL@ type for floats, where add would be defined as -(+.), for example. In our matrix functor, we will then use the functions defined -in the signature below (for example, add instead of (+.)). Also, we have -included some generating functions to help us with testing. - -{\annot{For the updated view of this, please take a look at EltsI.ml}} -\begin{verbatim} -module type ORDERED_AND_OPERATIONAL -sig - type t - - val zero : t - val one: t - val compare : t -> t -> order - - (* Converts a t to a string *) - val to_string : t -> string - - val add: t -> t -> t - val subtract: t -> t -> t - val multiply: t -> t -> t - val divide: t -> t -> t - - (* For testing *) - (* Prints a t *) - val print: t -> unit - - (* Generates the same t each time when called *) - val generate: unit -> t - - (* Generates a t greater than the argument passed in *) - val generate_gt: t -> unit -> t - - (* Generates a t less than the argument passed in *) - val generate_lt: t -> unit -> t - - (* Generates a t in between the two arguments. Returns none if none exist *) - val generate_between: t -> t -> unit -> t option -end -\end{verbatim} - -The following signature will provide all the necessary functions that can be -performed on a matrix. We include basic matrix operations like additions, -multiplications, scalar multiplication, and then standard operations like -finding the inverse of a matrix, finding its trace, determinant, eigenvalues, -transpose, norm, and how to row reduce the matrix. - -Below each function, in comments, is a brief description of how we will -implement the function. - -{\annot{For the updated view of this, please take a look at MatrixI.ml}} -\begin{verbatim} -module type MATRIX = -sig - exception NonSquare - - type elt - type matrix - (* Type of this is unknown, but will probably be represented using Ocaml’s - * built-in Arrays *) - (* empty matrix *) - val empty - - (* Takes a list of lists and converts that to a matrix *) - val from_list : (elt list list) -> matrix - (* Will implement using nested match statements *) - - (* Scales every element in the matrix by another elt *) - val scale : matrix -> elt -> matrix - (* Will implement by iterating through the matrix and scaling each element *) - - (* Adds two matrices. They must have the same dimensions *) - val add : matrix -> matrix -> matrix - (* Will add the elements elementwise and construct a new matrix *) - - (* Multiplies two matrices. If the matrices have dimensions m x n and p x q, n - * and p must be equal, and the resulting matrix will have dimension m x q *) - val mult: matrix -> matrix -> matrix - (* Will take the dot product of the nth row of the first matrix and the jth - * column of the second matrix to create the n,j th entry of the resultant *) - - (* Returns the row reduced form of a matrix *) - val row_reduce: matrix -> matrix - (* We will implement the algorithm found in the link above *) - - (* Returns the inverse of a matrix *) - val inverse: matrix -> matrix - (* Will implement this based on the specification in the Algorithms book *) - - (* Returns the norm of the matrix *) - val norm: matrix -> elt - - (* Transposes a matrix. If the input has dimensions m x n, the output will - * have dimensions n x m *) - val transpose: matrix -> matrix - (* Will basically ``flip'' the indices of the input matrix - - (* Returns the trace of the matrix *) - val trace: matrix -> elt - (* Will check if the matrix is square, then sum up all the elements along its - * diagonal *) - - (* Returns the determinant of the matrix *) - val det: matrix -> elt - (* Will implement this algorithm based on a description in Hubbard. Involves - * column reducing the input (or row-reducing the transpose) and then keeping - * track of the operations to build a sequence of coefficients to multiply *) - - (* Returns a list of eigenvalues and eigenvectors of a matrix *) - val eigen: matrix -> (elt *matrix) list option - (* Calculates successive powers of the input matrix, each multiplied by the - * same basis vector. Generates a polynomial and solves for zeros, which - * yields eigenvalues. Repeat for all basis vectors *) - - (* Takes a string and builds a matrix from it *) - from_string : string -> matrix - (* We will have some way to express matrices using strings, and then we will - * parse the string to give the matrix *) - - (* Prints out the contents of a matrix *) - print :matrix -> unit - (* Iterate through the matrix and print each element *) -end -\end{verbatim} - -The exception exists in case the user tries to perform operations, such as trace -or determinant, which require a square matrix. Matrices have their own type, and -they contain elements of the type \verb@elt@. - -We will probably need a helper function to raise matrices to powers and to solve -for zeros of a polynomial if we want to implement the \verb@eigen@ function. - -\subsection{Simplex Algorithm} - -The simplex solves linear programs, which are basically linear inequalities. -Given inequalities, using the simplex algorithm, we can find the maximum value -of a given expression. What the simplex algorithm does is that it takes in -several linear inequalities and writes them so that they are equalities. For -instance $5x + 7y \leq 10$ would be $5x + 7y + s = 10$ where $s \geq 0$ is the -slack variable. Then it takes all these equations and puts them a a matrix where -each entry represents the coefficient of the variable in the expression. The -first row represents the equation we are optimizing. The final column is the -constant values of each equation. Then we will do a series of computations so -that the top row has all nonnegative values. This can be done using the -functions from the Matrix module. - -The steps are as follows: - -Look along the first row and find the smallest negative element. Look at the -column of that entry and compute the value/element of each row. If any element -is zero, we can ignore it. Otherwise take the element with the smallest -value/element and make this the pivot point. Reduce the column so that the pivot -point in 1 and any other entry in the column is zero. Repeat the process again -until the top row consists of only nonnegative values. The first entry in the -last column is the maximum. In addition, for any column consisting of all zeroes -and a one, we can find the value of the corresponding variable by looking at the -value of the row containing the one. If the column doesn't consists of zeroes -and a one, then the value of the corresponding variable is zero. - -This algorithm also works with variables that are negative, inequalities where -the unknowns is greater than a constant, and minimum of the equation. This can -be easily done by either multiplying the whole inequality by -1 or by negating -the variable. If calculating the minimum or negating the variable, we can revert -them later. - -{\annot{For the updated view of this, please take a look at SimplexI.ml}} -\begin{verbatim} -module type simplex -sig - type elt - type linear_equation - type constraint - - val pivot : matrix -> matrix - val constraint_from_string : string -> constraint - val linear_from_string : string -> constraint -end -\end{verbatim} - -These are helper functions that we think will be used: - -\verb@no_neg : matrix -> bool@ - -Check the top row has no negative elements, which can be done by using the next -helper function - -\verb@min_of_row : row -> elt*(elt ref)@ - -A helper for the row reduction which will find the minimum element and its -location. -{\annot{We ended up needing many more helper functions than those shown here.}} - -\section{What's Next} - -We will read more into the algorithms and get a basic understanding of how they -all work so we can better judge the difficulty of each algorithm. We will also -read more into OCaml to see how to put together a project with multiple files -(since we have always relied on staff’s frameworks for each pset). We will -working through the CS50 appliance. We already have setup a git repository on -Github (\url{https://github.com/Fantastic-four}), and most of us are fairly -familiar with version control. - -\end{document} +\documentclass[letterpaper,11pt]{article} +% Change the header if you're going to change this. + +% Possible packages - uncomment to use +%\usepackage{amsmath} % Needed for math stuff +%\usepackage{lastpage} % Finds last page +%\usepackage{amssymb} % Needed for some math symbols +%\usepackage{graphicx} % Needed for graphics +\usepackage[usenames,dvipsnames]{xcolor} % Needed for graphics and color +%\usepackage{setspace} % Needed to set single/double spacing +\usepackage{hyperref} + +%Suggested by TeXworks +\usepackage[utf8]{inputenc} % set input encoding (not needed with XeLaTeX) +\usepackage{verbatim} % adds environment for commenting out blocks of text & for better verbatim + +% Sets the page margins to 1 inch each side +\usepackage[margin=1in]{geometry} + +% Uncomment this section if you wish to have a header. +\usepackage{fancyhdr} +\pagestyle{fancy} +\renewcommand{\headrulewidth}{0.5pt} % customise the layout... +\lhead{CS51 Draft Spec} \chead{} \rhead{May 5th, 2013} +\lfoot{} \cfoot{\thepage} \rfoot{} + +\frenchspacing +\setlength\parindent{0pt} +\setlength\parskip{2ex} + +\newcommand{\annot}[1]{\textbf{\color{BrickRed} [#1]}} + +\begin{document} +\title{CS 51 Final Project Functional Specification{\annot{Matrix Module and Simplex Algorithm}} +---{\annot{Annotated!}}} +\author{ +Luis Perez $|$ luisperez@college.harvard.edu \\ +Andy Shi $|$ andyshi@college.harvard.edu \\ +Zihao Wang $|$ zihaowang01@college.harvard.edu \\ +Ding Zhou $|$ dzhou@college.harvard.edu +} +\date{May 5, 2013} +\maketitle + +%Put your stuff here +{\annot{Note for the reader: This is the annotated version of our CS 51 +Final Project Functional Specifications. All annotations will appear in this color.}}. + +\section{Brief Overview} + +Our project seeks to implement an efficient matrix library for Ocaml. The end +goal is to utilize this matrix library in conjunction with the simplex algorithm +to find solution {\annot{Linear Programs}}. In order to accomplish this, +there will be multiple small checkpoints. We need to implement and efficient +matrix representation (so as to not use the naive implementation of lists of +lists){\annot{(Successfully accomplished. Our implementation uses array for efficienty)}}. + We also hope to provide a useful interface that allows the user to input +arbitrarily large matrices effectively (possible by mapping over the entries of +the matrix){\annot{The user can input matrices with a simple text file composed of comma +seperated values}}. One big challenge for the matrix library will be the implementation +of efficient algorithms for multiplication and row reduction. We expect the +multiplication to initially consists of the standard algorithm, but will +eventually grow into either block multiplication or Strassen {\annot{We were not +able to implement more sophisticated multiplication algorithms, but it seems that the current +algorithm works fairly well}}. We might even implement some way to chose between the +algorithms depending on the inputted matrices, but that would not an essential feature +{\annot{This feature can still easily be implemented, but there simply was no time. We +already knew this, though, by the time we wrote our Technical Specs}}. + +The row reduction algorithm will certainly be one of the most fundamental +algorithms in our matrix library implementation. From our current understanding +of simplex, it seems like it will heavily rely on row reduction {\annot{The +simplex algorithm ended up relying on a different type of row reduction, so we +couldn't just plug in a matrix into our row reduction algorithm. We had to do a +lot of extra work to get simplex. Though, many of the helper functions for row +reduction ended up being extremely useful for Simplex, so this was not as +pointless as some might think}}. The only algorithm we can come up with so far +is Gauss-Jordan, though we have found some algorithms which are more suited for +computation by computers with imprecise storage of floats {\annot{In fact, we +just had to be careful which element to choose as a pivot. Futhermore, we ended +up implementing arbitrary precision, via Ocaml's num library, so this became a +non-issue}}. + +The simplex part of the problem is the main goal, but if we can get down the +matrix module, then it appears as if the simplex algorithm should be +straightforward {\annot{This was definitely not the case}}. Eventually, we wish +to implement a simplex module that allows the input of arbitrary linear +equations with linear constraints, and then we parse this information into the +corresponding matrix and finally perform the simplex algorithm on that matrix in +order to return a solution to the specified problem.{\annot{This is fully functionalt. +Input can be made in the form of a text file with a linear program. The algorithm +runs and not only returns the right values for simple cases, but also checks for +cases that might initially appear to be unsolvable but aren't, and for cases which +are unbounded.}} + +Other tentative goals are to implement matrices that can accept either +arbitrarily large numbers (bignums) or arbitrarily precise floats +(\verb@int list * int stream@). It would also be interesting, though not +necessary, to implement multiple algorithms for each operation and choose from +them the one that is the best for a specific input data {\annot{We were able to +implement arbitrary precisions and arbitrarily large sizes with the OCaml Nums +library}}. + +\section{Feature List} +\begin{enumerate} + +\item Find a way to represent a matrix (arrays or list list or other +representation)---will likely use an array representation using the Ocaml array +module (\url{http://caml.inria.fr/pub/docs/manual-ocaml/libref/Array.html}) +{\annot{Implemented as specified}} + +\item Make a matrix from a list of lists (this is how the user will input the +data of the matrix). {\annot{Implemented as specified}} + +\item Scalar multiplication---as described here +(\url{http://www.purplemath.com/modules/mtrxmult.htm}) +{\annot{Implemented as specified}} + +\item Matrix addition---as described here +(\url{http://www.purplemath.com/modules/mtrxmult.htm}) +{\annot{Implemented as specified}} + +\item Standard Matrix multiplication +(\url{http://www.purplemath.com/modules/mtrxmult.htm}) +{\annot{Implemented as specified}} + +\item Row Reduction (Gauss-Jordan). We will use the following resource to help +us calculate row reduction when the entries of the matrix are stored with +limited precision +(\url{http://thejuniverse.org/PUBLIC/LinearAlgebra/LOLA/rowRed/var.html}) +{\annot{Implemented as specified}} + +\item Matrix inverse / Row reduction with LU Decomposition - Description can be +found here (\url{http://www.math.ust.hk/~macheng/math111/LU_Decomposition.pdf}) +{\annot{To find the inverse, we augmented the original matrix with an identity +matrix, then row-reduced the augmented matrix}} + +\item{\annot{Make a matrix from a text file. The user is now able to read in a +text file and create a matrix.}} + +\item{\annot{For Information of further functionality, take a look at MatrixI.ml.}} + +\item Simplex Algorithm (Described in Algorithms book, in addition to the +resources provided below) +{\annot{Implemented mostly as specified, with some modifications from the +Wikipedia article on simplex}} + +\item Matrix norm +{\annot{Did not implement since it was not necessary for Simplex.}} + +\item Matrix transpose +{\annot{Implemented as specified.}} + +\item Trace---This is a typical trace (sum of the entries on the diagonal) +{\annot{Implemented as specified.}} + +\item Determinant---found using algorithm described in Section 4.8 of Hubbard's +\emph{Vector Calculus, Linear Algebra, and Differential Forms: A Unified +Approach} (Math 23a/b textbook) +{\annot{Implemented by decomposing the matrix into upper and lower triangular +matrices, then multiplying the entries along the diagonals for both and +multiplying the two results}} + +\item Eigenvectors / Eigenvalues, calculated also by a method in Hubbard's book +(this might not be necessary to carry out the simplex algorithm and we might +abandon it if we can't figure out how to implement it in a timely manner) +{\annot{Did not implement since it was not necssary for Simplex.}} + +\end{enumerate} + +\section{Draft Technical Specification} + +\subsection{Matrices} + +Our matrix module will actually be a functor that takes in an +\verb@ORDERED_AND_OPERATIONAL@ module (from ps4 and moogle) so we can have +matrices of ints, floats, or floats as streams. +{\annot{This actually ended up being very helpful since at the end we we able +to use OCaml's Num Library in order to achieve arbitrary precision and size.}} + +The following type definition (a la ps4) will give us a way to talk about order +for different data types. + +\begin{verbatim} +type order = Equal | Less | Greater +\end{verbatim} + +The following signature will allow us to provide matrix functionality on a +variety of data types. For example, we could make an +\verb@ORDERED_AND_OPERATIONAL@ type for floats, where add would be defined as +(+.), for example. In our matrix functor, we will then use the functions defined +in the signature below (for example, add instead of (+.)). Also, we have +included some generating functions to help us with testing. + +{\annot{For the updated view of this, please take a look at EltsI.ml}} +\begin{verbatim} +module type ORDERED_AND_OPERATIONAL +sig + type t + + val zero : t + val one: t + val compare : t -> t -> order + + (* Converts a t to a string *) + val to_string : t -> string + + val add: t -> t -> t + val subtract: t -> t -> t + val multiply: t -> t -> t + val divide: t -> t -> t + + (* For testing *) + (* Prints a t *) + val print: t -> unit + + (* Generates the same t each time when called *) + val generate: unit -> t + + (* Generates a t greater than the argument passed in *) + val generate_gt: t -> unit -> t + + (* Generates a t less than the argument passed in *) + val generate_lt: t -> unit -> t + + (* Generates a t in between the two arguments. Returns none if none exist *) + val generate_between: t -> t -> unit -> t option +end +\end{verbatim} + +The following signature will provide all the necessary functions that can be +performed on a matrix. We include basic matrix operations like additions, +multiplications, scalar multiplication, and then standard operations like +finding the inverse of a matrix, finding its trace, determinant, eigenvalues, +transpose, norm, and how to row reduce the matrix. + +Below each function, in comments, is a brief description of how we will +implement the function. + +{\annot{For the updated view of this, please take a look at MatrixI.ml}} +\begin{verbatim} +module type MATRIX = +sig + exception NonSquare + + type elt + type matrix + (* Type of this is unknown, but will probably be represented using Ocaml’s + * built-in Arrays *) + (* empty matrix *) + val empty + + (* Takes a list of lists and converts that to a matrix *) + val from_list : (elt list list) -> matrix + (* Will implement using nested match statements *) + + (* Scales every element in the matrix by another elt *) + val scale : matrix -> elt -> matrix + (* Will implement by iterating through the matrix and scaling each element *) + + (* Adds two matrices. They must have the same dimensions *) + val add : matrix -> matrix -> matrix + (* Will add the elements elementwise and construct a new matrix *) + + (* Multiplies two matrices. If the matrices have dimensions m x n and p x q, n + * and p must be equal, and the resulting matrix will have dimension m x q *) + val mult: matrix -> matrix -> matrix + (* Will take the dot product of the nth row of the first matrix and the jth + * column of the second matrix to create the n,j th entry of the resultant *) + + (* Returns the row reduced form of a matrix *) + val row_reduce: matrix -> matrix + (* We will implement the algorithm found in the link above *) + + (* Returns the inverse of a matrix *) + val inverse: matrix -> matrix + (* Will implement this based on the specification in the Algorithms book *) + + (* Returns the norm of the matrix *) + val norm: matrix -> elt + + (* Transposes a matrix. If the input has dimensions m x n, the output will + * have dimensions n x m *) + val transpose: matrix -> matrix + (* Will basically ``flip'' the indices of the input matrix + + (* Returns the trace of the matrix *) + val trace: matrix -> elt + (* Will check if the matrix is square, then sum up all the elements along its + * diagonal *) + + (* Returns the determinant of the matrix *) + val det: matrix -> elt + (* Will implement this algorithm based on a description in Hubbard. Involves + * column reducing the input (or row-reducing the transpose) and then keeping + * track of the operations to build a sequence of coefficients to multiply *) + + (* Returns a list of eigenvalues and eigenvectors of a matrix *) + val eigen: matrix -> (elt *matrix) list option + (* Calculates successive powers of the input matrix, each multiplied by the + * same basis vector. Generates a polynomial and solves for zeros, which + * yields eigenvalues. Repeat for all basis vectors *) + + (* Takes a string and builds a matrix from it *) + from_string : string -> matrix + (* We will have some way to express matrices using strings, and then we will + * parse the string to give the matrix *) + + (* Prints out the contents of a matrix *) + print :matrix -> unit + (* Iterate through the matrix and print each element *) +end +\end{verbatim} + +The exception exists in case the user tries to perform operations, such as trace +or determinant, which require a square matrix. Matrices have their own type, and +they contain elements of the type \verb@elt@. + +We will probably need a helper function to raise matrices to powers and to solve +for zeros of a polynomial if we want to implement the \verb@eigen@ function. + +\subsection{Simplex Algorithm} + +The simplex solves linear programs, which are basically linear inequalities. +Given inequalities, using the simplex algorithm, we can find the maximum value +of a given expression. What the simplex algorithm does is that it takes in +several linear inequalities and writes them so that they are equalities. For +instance $5x + 7y \leq 10$ would be $5x + 7y + s = 10$ where $s \geq 0$ is the +slack variable. Then it takes all these equations and puts them a a matrix where +each entry represents the coefficient of the variable in the expression. The +first row represents the equation we are optimizing. The final column is the +constant values of each equation. Then we will do a series of computations so +that the top row has all nonnegative values. This can be done using the +functions from the Matrix module. + +The steps are as follows: + +Look along the first row and find the smallest negative element. Look at the +column of that entry and compute the value/element of each row. If any element +is zero, we can ignore it. Otherwise take the element with the smallest +value/element and make this the pivot point. Reduce the column so that the pivot +point in 1 and any other entry in the column is zero. Repeat the process again +until the top row consists of only nonnegative values. The first entry in the +last column is the maximum. In addition, for any column consisting of all zeroes +and a one, we can find the value of the corresponding variable by looking at the +value of the row containing the one. If the column doesn't consists of zeroes +and a one, then the value of the corresponding variable is zero. + +This algorithm also works with variables that are negative, inequalities where +the unknowns is greater than a constant, and minimum of the equation. This can +be easily done by either multiplying the whole inequality by -1 or by negating +the variable. If calculating the minimum or negating the variable, we can revert +them later. + +{\annot{For the updated view of this, please take a look at SimplexI.ml}} +\begin{verbatim} +module type simplex +sig + type elt + type linear_equation + type constraint + + val pivot : matrix -> matrix + val constraint_from_string : string -> constraint + val linear_from_string : string -> constraint +end +\end{verbatim} + +These are helper functions that we think will be used: + +\verb@no_neg : matrix -> bool@ + +Check the top row has no negative elements, which can be done by using the next +helper function + +\verb@min_of_row : row -> elt*(elt ref)@ + +A helper for the row reduction which will find the minimum element and its +location. +{\annot{We ended up needing many more helper functions than those shown here.}} + +\section{What's Next} + +We will read more into the algorithms and get a basic understanding of how they +all work so we can better judge the difficulty of each algorithm. We will also +read more into OCaml to see how to put together a project with multiple files +(since we have always relied on staff’s frameworks for each pset). We will +working through the CS50 appliance. We already have setup a git repository on +Github (\url{https://github.com/Fantastic-four}), and most of us are fairly +familiar with version control. + +\end{document} diff --git a/Elts.ml b/Elts.ml index c339e24..3c69197 100644 --- a/Elts.ml +++ b/Elts.ml @@ -1,159 +1,159 @@ -(* At the bottom, setting the Elts module equal to one of the following uses - * that module for the creation of matrices and simplex *) - -module Floats : EltsI.ORDERED_AND_OPERATIONAL = -struct - - exception NonElt - - type t = float - - (* Need an epsilon to compare floats. We define the precion here *) - let epsilon = 0.000001 - - let zero = 0. - - let one = 1. - - (* We compare based on epsilon. We first find the distance - * between the two numbers and divide my the maximum to see - * if that's below epsilon. *) - let compare a b = - let a', b' = abs_float a, abs_float b in - if a' < epsilon && b' < epsilon then - Order.Equal - else - let diff = (a -. b) /. (max a' b') in - if abs_float diff < epsilon then Order.Equal - else if a < b then Order.Less - else if a > b then Order.Greater - else - raise (Failure "Error in compare.") - - let to_string = string_of_float - - let from_string (x: string) = - try - float_of_string x - with - | Failure ("float_of_string") -> raise NonElt - - let add = (+.) - - let subtract = (-.) - - let multiply = ( *. ) - - let divide = (/.) - - let print a = print_string ((to_string a) ^ "\n") - - let trim = int_of_float - - let generate () = 3. - - let generate_gt a () = a +. 1. - - let generate_lt a () = a -. 1. - - let generate_between a b () = if a > b then None else Some((a +. b)/. 2.) - - let generate_x x () = (x:t) - - let generate_random bound () = - let x = Random.float bound in - x -. mod_float x epsilon - - (************************ TESTS ********************************) - - let rec test_compare (times:int) : unit = - let random t = float_of_int (Random.int t - Random.int t) in - if times = 0 then () - else - let x, y = random times, random times in - match compare x y with - | Order.Equal -> assert(x = y) - | Order.Greater -> assert(x > y) - | Order.Less -> assert(x < y) - - - (* Honestly, probably don't need to test the other functions *) - let run_tests (times: int) : unit = - test_compare times ; - () -end - -module Num_floats : EltsI.ORDERED_AND_OPERATIONAL = -struct - - exception NonElt - - type t = Num.num - - let zero = Num.num_of_int 0 - - let one = Num.num_of_int 1 - - let compare a b = - if Num.compare_num a b = -1 then Order.Less - else if Num.compare_num a b = 0 then Order.Equal - else if Num.compare_num a b = 1 then Order.Greater - else - raise (Failure "Error in compare.") - - let to_string = Num.string_of_num - - let from_string (x: string) = - try - Num.num_of_string x - with - | Failure ("num_of_string") -> raise NonElt - - let add = Num.add_num - - let subtract = Num.sub_num - - let multiply = Num.mult_num - - let divide = Num.div_num - - let generate () = Num.num_of_int 3 - - let generate_gt a () = Num.add_num a one - - let generate_lt a () = Num.sub_num a (subtract zero one) - - let generate_between a b () = - if Num.gt_num a b then None - else Some (Num.div_num (Num.add_num a b) (Num.num_of_int 2)) - - let generate_x (x:float) () = - let new_x = int_of_float x in Num.num_of_int new_x - - (* This actually only generates integers. *) - let generate_random (bound: float) () = - let new_bound = int_of_float bound in - Num.num_of_int (Random.int new_bound) - - let print a = print_string ((to_string a) ^ "\n") - - let rec test_compare (times:int) : unit = - let random t = Num.num_of_int (Random.int t - Random.int t) in - if times = 0 then () - else - let x, y = random times, random times in - match compare x y with - | Order.Equal -> assert(x = y) - | Order.Greater -> assert(x > y) - | Order.Less -> assert(x < y) - - - let run_tests (times: int) : unit = - test_compare times ; - () -end - - -(*** HERE IS WHERE YOU PICK THE MODULE *) -(*module Elts = Floats *) -module Elts = Num_floats +(* At the bottom, setting the Elts module equal to one of the following uses + * that module for the creation of matrices and simplex *) + +module Floats : EltsI.ORDERED_AND_OPERATIONAL = +struct + + exception NonElt + + type t = float + + (* Need an epsilon to compare floats. We define the precion here *) + let epsilon = 0.000001 + + let zero = 0. + + let one = 1. + + (* We compare based on epsilon. We first find the distance + * between the two numbers and divide my the maximum to see + * if that's below epsilon. *) + let compare a b = + let a', b' = abs_float a, abs_float b in + if a' < epsilon && b' < epsilon then + Order.Equal + else + let diff = (a -. b) /. (max a' b') in + if abs_float diff < epsilon then Order.Equal + else if a < b then Order.Less + else if a > b then Order.Greater + else + raise (Failure "Error in compare.") + + let to_string = string_of_float + + let from_string (x: string) = + try + float_of_string x + with + | Failure ("float_of_string") -> raise NonElt + + let add = (+.) + + let subtract = (-.) + + let multiply = ( *. ) + + let divide = (/.) + + let print a = print_string ((to_string a) ^ "\n") + + let trim = int_of_float + + let generate () = 3. + + let generate_gt a () = a +. 1. + + let generate_lt a () = a -. 1. + + let generate_between a b () = if a > b then None else Some((a +. b)/. 2.) + + let generate_x x () = (x:t) + + let generate_random bound () = + let x = Random.float bound in + x -. mod_float x epsilon + + (************************ TESTS ********************************) + + let rec test_compare (times:int) : unit = + let random t = float_of_int (Random.int t - Random.int t) in + if times = 0 then () + else + let x, y = random times, random times in + match compare x y with + | Order.Equal -> assert(x = y) + | Order.Greater -> assert(x > y) + | Order.Less -> assert(x < y) + + + (* Honestly, probably don't need to test the other functions *) + let run_tests (times: int) : unit = + test_compare times ; + () +end + +module Num_floats : EltsI.ORDERED_AND_OPERATIONAL = +struct + + exception NonElt + + type t = Num.num + + let zero = Num.num_of_int 0 + + let one = Num.num_of_int 1 + + let compare a b = + if Num.compare_num a b = -1 then Order.Less + else if Num.compare_num a b = 0 then Order.Equal + else if Num.compare_num a b = 1 then Order.Greater + else + raise (Failure "Error in compare.") + + let to_string = Num.string_of_num + + let from_string (x: string) = + try + Num.num_of_string x + with + | Failure ("num_of_string") -> raise NonElt + + let add = Num.add_num + + let subtract = Num.sub_num + + let multiply = Num.mult_num + + let divide = Num.div_num + + let generate () = Num.num_of_int 3 + + let generate_gt a () = Num.add_num a one + + let generate_lt a () = Num.sub_num a (subtract zero one) + + let generate_between a b () = + if Num.gt_num a b then None + else Some (Num.div_num (Num.add_num a b) (Num.num_of_int 2)) + + let generate_x (x:float) () = + let new_x = int_of_float x in Num.num_of_int new_x + + (* This actually only generates integers. *) + let generate_random (bound: float) () = + let new_bound = int_of_float bound in + Num.num_of_int (Random.int new_bound) + + let print a = print_string ((to_string a) ^ "\n") + + let rec test_compare (times:int) : unit = + let random t = Num.num_of_int (Random.int t - Random.int t) in + if times = 0 then () + else + let x, y = random times, random times in + match compare x y with + | Order.Equal -> assert(x = y) + | Order.Greater -> assert(x > y) + | Order.Less -> assert(x < y) + + + let run_tests (times: int) : unit = + test_compare times ; + () +end + + +(*** HERE IS WHERE YOU PICK THE MODULE *) +(*module Elts = Floats *) +module Elts = Num_floats diff --git a/EltsI.ml b/EltsI.ml index 5724172..8c475fd 100644 --- a/EltsI.ml +++ b/EltsI.ml @@ -1,60 +1,60 @@ -module type ORDERED_AND_OPERATIONAL = -sig - - (* Exception for from_string. Is raised when from_string is passed something - * that is not an elt *) - exception NonElt - - type t - - (* The zero element *) - val zero : t - - (* The one element *) - val one: t - - (* ts must be comparable *) - val compare : t -> t -> Order.order - - (* Converts a t to a string *) - val to_string : t -> string - - (* Converts a string to a t *) - val from_string : string -> t - - (* Basic mathematical operations must be possible *) - val add: t -> t -> t - - val subtract: t -> t -> t - - val multiply: t -> t -> t - - val divide: t -> t -> t - - (* Prints a t *) - val print: t -> unit - - (************* For testing ***************) - (* Generates the same t each time when called *) - val generate: unit -> t - - (* Generates a t greater than the argument passed in *) - val generate_gt: t -> unit -> t - - (* Generates a t less than the argument passed in *) - val generate_lt: t -> unit -> t - - (* Generates a t in between the two arguments. Returns none if none exist *) - val generate_between: t -> t -> unit -> t option - - (* Special test function specifically for float. When passed in a float x, - * generates an abstract float *) - val generate_x: float -> unit -> t - - (* Special function for testing. Generates a random t between the zero element - * and the bound (the argument) inclusive *) - val generate_random : float -> unit -> t - - val run_tests : int -> unit - -end +module type ORDERED_AND_OPERATIONAL = +sig + + (* Exception for from_string. Is raised when from_string is passed something + * that is not an elt *) + exception NonElt + + type t + + (* The zero element *) + val zero : t + + (* The one element *) + val one: t + + (* ts must be comparable *) + val compare : t -> t -> Order.order + + (* Converts a t to a string *) + val to_string : t -> string + + (* Converts a string to a t *) + val from_string : string -> t + + (* Basic mathematical operations must be possible *) + val add: t -> t -> t + + val subtract: t -> t -> t + + val multiply: t -> t -> t + + val divide: t -> t -> t + + (* Prints a t *) + val print: t -> unit + + (************* For testing ***************) + (* Generates the same t each time when called *) + val generate: unit -> t + + (* Generates a t greater than the argument passed in *) + val generate_gt: t -> unit -> t + + (* Generates a t less than the argument passed in *) + val generate_lt: t -> unit -> t + + (* Generates a t in between the two arguments. Returns none if none exist *) + val generate_between: t -> t -> unit -> t option + + (* Special test function specifically for float. When passed in a float x, + * generates an abstract float *) + val generate_x: float -> unit -> t + + (* Special function for testing. Generates a random t between the zero element + * and the bound (the argument) inclusive *) + val generate_random : float -> unit -> t + + val run_tests : int -> unit + +end diff --git a/Helpers.ml b/Helpers.ml index 007dfdb..6980052 100644 --- a/Helpers.ml +++ b/Helpers.ml @@ -1,16 +1,16 @@ -(* Takes in a string and a separator, and separates the string into a list of - * substrings where each substring is between two separators or between a - * separator and the beginning/end of the string *) -let explode (s: string) (space: string) : string list = - let rec build (curr: string) (buffer: string) (lst: string list) : string list = - let len = String.length curr in - if len = 0 then buffer::lst - else - let c = String.sub curr (len - 1) 1 in - if len = 1 then (c ^ buffer)::lst - else - let s' = String.sub curr 0 (len - 1) in - if c = space then build s' "" (buffer::lst) - else build s' (c ^ buffer) lst in - build (String.trim s) "" [] - +(* Takes in a string and a separator, and separates the string into a list of + * substrings where each substring is between two separators or between a + * separator and the beginning/end of the string *) +let explode (s: string) (space: string) : string list = + let rec build (curr: string) (buffer: string) (lst: string list) : string list = + let len = String.length curr in + if len = 0 then buffer::lst + else + let c = String.sub curr (len - 1) 1 in + if len = 1 then (c ^ buffer)::lst + else + let s' = String.sub curr 0 (len - 1) in + if c = space then build s' "" (buffer::lst) + else build s' (c ^ buffer) lst in + build (String.trim s) "" [] + diff --git a/Interface.ml b/Interface.ml index d520d51..7da3617 100644 --- a/Interface.ml +++ b/Interface.ml @@ -1,52 +1,52 @@ -open Elts -open Matrix -open MatrixI -open SimplexI - -(* Testing the Matrix Library *) -let test times = - let _ = EltMatrix.run_tests times in - let _ = Elts.run_tests times in - let _ = Simplex.run_tests times in - print_string "\nIf no errors above, then all tests passed!\n" - -(* Prints a solution *) -let print_solution (e,p) : unit = - let _ = print_string "Solved!\n\nOptimal value: " in - let _ = Elts.print e in - let _ = print_string "\n\nAchieved at the point: " in - let _ = Simplex.print_point p in - print_string "\n\n" - - - -(* Parse command-line arguments. Parses filename *) -let parse_args () : unit = - let usage s = - let message = "Possible commands:\n'simplex' {file}\n'run tests' {times}\n" - in - Printf.printf ("usage: %s arguments.\nError-> %s\n") - Sys.argv.(0) s; print_string message; exit 1 in - if Array.length Sys.argv <> 3 then usage "Incorrect number of arguments."; - match String.lowercase (Sys.argv.(1)) with - | "run tests" | "run_tests" -> - (try - let times = int_of_string Sys.argv.(2) in - test times - with - | Failure s -> usage s) - | "solve"-> - (try - let filename = Sys.argv.(2) in - (match Simplex.load_file filename with - | None -> (print_string "\nThis system has no feasable solution.\n") - | Some sys -> - let _ = print_string "\nSolving your system....\n\n" in - match Simplex.solve sys with - | None -> print_string "This system is unbounded. You can increase/decrease it as you please!\n" - | Some solution -> print_solution solution) - with - | Sys_error e -> usage e) - | _ -> usage "Incorrect inputs." -;; - +open Elts +open Matrix +open MatrixI +open SimplexI + +(* Testing the Matrix Library *) +let test times = + let _ = EltMatrix.run_tests times in + let _ = Elts.run_tests times in + let _ = Simplex.run_tests times in + print_string "\nIf no errors above, then all tests passed!\n" + +(* Prints a solution *) +let print_solution (e,p) : unit = + let _ = print_string "Solved!\n\nOptimal value: " in + let _ = Elts.print e in + let _ = print_string "\n\nAchieved at the point: " in + let _ = Simplex.print_point p in + print_string "\n\n" + + + +(* Parse command-line arguments. Parses filename *) +let parse_args () : unit = + let usage s = + let message = "Possible commands:\n'simplex' {file}\n'run tests' {times}\n" + in + Printf.printf ("usage: %s arguments.\nError-> %s\n") + Sys.argv.(0) s; print_string message; exit 1 in + if Array.length Sys.argv <> 3 then usage "Incorrect number of arguments."; + match String.lowercase (Sys.argv.(1)) with + | "run tests" | "run_tests" -> + (try + let times = int_of_string Sys.argv.(2) in + test times + with + | Failure s -> usage s) + | "solve"-> + (try + let filename = Sys.argv.(2) in + (match Simplex.load_file filename with + | None -> (print_string "\nThis system has no feasable solution.\n") + | Some sys -> + let _ = print_string "\nSolving your system....\n\n" in + match Simplex.solve sys with + | None -> print_string "This system is unbounded. You can increase/decrease it as you please!\n" + | Some solution -> print_solution solution) + with + | Sys_error e -> usage e) + | _ -> usage "Incorrect inputs." +;; + diff --git a/Main.ml b/Main.ml index cea3c2d..a1b8076 100644 --- a/Main.ml +++ b/Main.ml @@ -1,25 +1,25 @@ -(* PROJECT: OCaml Matrix Library/Simplex - * - * NAMES: - * - * Partner 1's name: Andy Shi - * Partner 1's github username: shiandy - * - * Partner 2's name: Ding Zhou - * Partner 2's github username: dzhou94 - * - * Partner 3's name: Jason Wang - * Partner 3's github username: jasonwang99 - * - * Partner 4's name: Luis Perez - * Partner 4's github username: kandluis - * - * Github Organization name: Fantastic-four - * - * Github Project URL: github.com/Fantastic-four/ocaml-matrix - * SEAS Project URL: code.seas.harvard.edu/ocaml-matrix/ocaml-matrix - * Video URL: - *) - -(* parses the arguments from the command-line. Interface loads everything up *) +(* PROJECT: OCaml Matrix Library/Simplex + * + * NAMES: + * + * Partner 1's name: Andy Shi + * Partner 1's github username: shiandy + * + * Partner 2's name: Ding Zhou + * Partner 2's github username: dzhou94 + * + * Partner 3's name: Jason Wang + * Partner 3's github username: jasonwang99 + * + * Partner 4's name: Luis Perez + * Partner 4's github username: kandluis + * + * Github Organization name: Fantastic-four + * + * Github Project URL: github.com/Fantastic-four/ocaml-matrix + * SEAS Project URL: code.seas.harvard.edu/ocaml-matrix/ocaml-matrix + * Video URL: + *) + +(* parses the arguments from the command-line. Interface loads everything up *) Interface.parse_args () ;; \ No newline at end of file diff --git a/Makefile b/Makefile index 2bb0ea4..13759f9 100644 --- a/Makefile +++ b/Makefile @@ -1,43 +1,43 @@ -# This program - -PROG = simplex - -# Setup - -LIBS = \ - nums.cma - -CAMLC = ocamlc -CAMLDOC = ocamldoc -CAMLFLAGS = -g - -%.cmo: %.ml - $(CAMLC) $(CAMLFLAGS) -c $< - -# Source and Object files -SOURCES = \ - Order.ml EltsI.ml Elts.ml \ - Helpers.ml MatrixI.ml Matrix.ml \ - SimplexI.ml Interface.ml Main.ml - -OBJECTS = $(SOURCES:.ml=.cmo) - -# Final Program - -$(PROG): $(OBJECTS) - $(CAMLC) $(CAMLFLAGS) $(LIBS) $(OBJECTS) -o $(PROG) - -# DocGen - -doc: $(OBJECTS) - $(CAMLDOC) -html $(SOURCES) - -# Other - -all: $(PROG) - -clean: - rm -rf *.cmo *.cmi *.html *.css $(PROG) - -.DEFAULT_GOAL := $(PROG) -.PHONY: doc build run clean +# This program + +PROG = simplex + +# Setup + +LIBS = \ + nums.cma + +CAMLC = ocamlc +CAMLDOC = ocamldoc +CAMLFLAGS = -g + +%.cmo: %.ml + $(CAMLC) $(CAMLFLAGS) -c $< + +# Source and Object files +SOURCES = \ + Order.ml EltsI.ml Elts.ml \ + Helpers.ml MatrixI.ml Matrix.ml \ + SimplexI.ml Interface.ml Main.ml + +OBJECTS = $(SOURCES:.ml=.cmo) + +# Final Program + +$(PROG): $(OBJECTS) + $(CAMLC) $(CAMLFLAGS) $(LIBS) $(OBJECTS) -o $(PROG) + +# DocGen + +doc: $(OBJECTS) + $(CAMLDOC) -html $(SOURCES) + +# Other + +all: $(PROG) + +clean: + rm -rf *.cmo *.cmi *.html *.css $(PROG) + +.DEFAULT_GOAL := $(PROG) +.PHONY: doc build run clean diff --git a/Matrix.ml b/Matrix.ml index eae1e39..f7f884b 100644 --- a/Matrix.ml +++ b/Matrix.ml @@ -1,802 +1,802 @@ -open Order -open Elts - -(* Functor so we can Abstract away! *) -module MakeMatrix (C: EltsI.ORDERED_AND_OPERATIONAL) : - (MatrixI.MATRIX with type elt = C.t) = -struct - - (*************** Exceptions ***************) - - exception NonSquare - exception ImproperDimensions - - (*************** End Exceptions ***************) - - (*************** Types ***************) - - type elt = C.t - - (* A matrix is a pair of dimension (n x p) and a array of arrays - * the first array is the row (n) and the second the column (p) *) - type matrix = (int * int) * (elt array array) - - (*************** End Types ***************) - - (*************** Base Functions ***************) - - (* catching negative dimensions AND 0 dimensions and too large - * of a dimension so we don't have to worry about it later *) - let empty (rows: int) (columns: int) : matrix = - if rows > 0 && columns > 0 then - try - let m = Array.make_matrix rows columns C.zero in ((rows,columns),m) - with Invalid_argument "index out of bounds" -> - raise ImproperDimensions - else (* dimension is negative or 0 *) - raise ImproperDimensions - - (*************** End Base Functions ***************) - - (*************** Helper Functions ***************) - - (* get's the nth row of a matrix and returns (r, row) where r is the length - * of the row and row is a COPY of the original row. For example, calling - * calling get_row m 1 will return (3, |1 3 4 |) - * ________ - * m = | 1 3 4 | - * |*2 5 6 | - *) - (* aside: we don't check whether n < 1 because of our matrix invariant *) - let get_row (((n,p),m): matrix) (row: int) : int * elt array = - if row <= n then - let row' = Array.map (fun x -> x) m.(row - 1) in - (p, row') - else - raise (Failure "Row out of bounds.") - - (* similar to get_row. For m, get_column m 1 will return (2, |1 2|) *) - let get_column (((n,p),m): matrix) (column: int) : int * elt array = - if column <= p then - begin - let column' = Array.make n C.zero in - for i = 0 to n - 1 do - column'.(i) <- m.(i).(column - 1) - done; - (n, column') - end - else - raise (Failure "Column out of bounds.") - - (* sets the nth row of the matrix m to the specified array a. - * This is done IN-PLACE. Therefore the function returns unit. You should - * nonetheless enfore immutability whenever possible. For a clarification on - * what nth row means, look at comment for get_row above. *) - let set_row (((n,p),m): matrix) (row: int) (a: elt array) : unit = - if row <= n then - begin - assert(Array.length a = p); - for i = 0 to p - 1 do - m.(row - 1).(i) <- a.(i) - done; - end - else - raise (Failure "Row out of bounds.") - - (* Similar to set_row but sets the nth column instead *) - let set_column (((n,p),m): matrix) (column: int) (a: elt array) : unit = - if column <= p then - begin - assert(Array.length a = n); - for i = 0 to n - 1 do - m.(i).(column - 1) <- a.(i) - done; - end - else - raise (Failure "Column out of bounds.") - - (* returns the ij-th element of a matrix (not-zero indexed) *) - let get_elt (((n,p),m): matrix) ((i,j): int*int) : elt = - if i <= n && j <= p then - m.(i - 1).(j - 1) - else - raise ImproperDimensions - - (* Changes the i,j-th element of a matrix to e. Is not zero-indexed, and - * changes the matrix in place *) - let set_elt (((n,p),m): matrix) ((i,j): int*int) (e: elt) : unit = - if i <= n && j <= p then - m.(i - 1).(j - 1) <- e - else - raise ImproperDimensions - - (* similar to map, but applies to function to the entire matrix - * Returns a new matrix *) - let map (f: elt -> elt) (mat: matrix) : matrix = - let (dim,m) = mat in - (dim, Array.map (Array.map f) m) - - (* Just some wrapping of Array.iter made for Matrices! *) - let iter (f: elt -> unit) (mat: matrix) : unit = - let dim,m = mat in - Array.iter (Array.iter f) m - - (* Just some wrapping of Array.iteri. Useful for pretty - * printing matrix. The index is (i,j). NOT zero-indexed *) - let iteri (f: int -> int -> elt -> unit) (mat: matrix) : unit = - let dim,m = mat in - Array.iteri (fun i row -> Array.iteri (fun j e -> f (i+1) (j+1) e) row) m - - (* folds over each row using base case u and function f *) - (* could be a bit more efficient? *) - let reduce (f: 'a -> elt -> 'a) (u: 'a) (((p,q),m): matrix) : 'a = - let total = ref u in - for i = 0 to p - 1 do - for j = 0 to q - 1 do - total := f (!total) m.(i).(j) - done; - done; - !total - - (* given two arrays, this will calculate their dot product *) - (* It seems a little funky, but this is done for efficiency's sake. - * In short, it tries to multiply each element by it's respective - * element until the one array is indexed out of bounds. If the - * other array is also out of bounds, then it returns their value. - * Otherwise, the arrays were the wrong size and raises ImproperDimension - - THE ABOVE COMMENT HAS NOT BEEN IMPLEMENTED - - Instead we calculate the length before starting - *) - let dot (v1: elt array) (v2: elt array) : elt = - let rec dotting (i: int) (total: elt) : elt = - if i = 0 then total - else - let curr = C.multiply v1.(i-1) v2.(i-1) in - dotting (i - 1) (C.add curr total) in - let len1, len2 = Array.length v1, Array.length v2 in - if len1 = len2 then dotting len1 C.zero - else raise ImproperDimensions - - (* function to expose the dimensions of a matrix *) - let get_dimensions (m: matrix) : (int * int) = - let ((x,y),m') = m in (x,y) - - (*************** End Helper Functions ***************) - - - (*************** Primary Matrix Functions ***************) - - (* scales a matrix by the appropriate factor *) - let scale (m: matrix) (sc: elt) : matrix = map (C.multiply sc) m - - (* Generates a matrix from a list of lists. The inners lists are the rows *) - let from_list (lsts : elt list list) : matrix = - let rec check_length (length: int) (lst: elt list) : int = - if List.length lst = length then length - else raise ImproperDimensions in - let p = List.length lsts in - match lsts with - | [] -> raise ImproperDimensions - | hd::tl -> - let len = List.length hd in - if List.fold_left check_length len tl = len then - ((p,len),Array.map Array.of_list (Array.of_list lsts)) - else - raise ImproperDimensions - - (* Adds two matrices. They must have the same dimensions *) - let add ((dim1,m1): matrix) ((dim2,m2): matrix) : matrix = - if dim1 = dim2 then - let n, p = dim1 in - let (dim', sum_m) = empty n p in - for i = 0 to n - 1 do - for j = 0 to p - 1 do - sum_m.(i).(j) <- C.add m1.(i).(j) m2.(i).(j) - done; - done; - (dim',sum_m) - else - raise ImproperDimensions - - - (* Multiplies two matrices. If the matrices have dimensions m x n and p x q, n - * and p must be equal, and the resulting matrix will have dimension n x q *) - let mult (matrix1: matrix) (matrix2: matrix) : matrix = - let ((m,n),m1), ((p,q),m2) = matrix1, matrix2 in - if n = p then - let (dim, result) = empty m q in - for i = 0 to m - 1 do - for j = 0 to q - 1 do - let (_,row), (_,column) = get_row matrix1 (i + 1), - get_column matrix2 (j + 1) in - result.(i).(j) <- dot row column - done; - done; - (dim,result) - else - raise ImproperDimensions - - (*************** Helper Functions for Row Reduce ***************) - - (* returns the index of the first non-zero elt in an array*) - let zero (arr: elt array) : int option = - let index = ref 1 in - let empty (i: int option) (e: elt) : int option = - match i, C.compare e C.zero with - | None, Equal -> (index := !index + 1; None) - | None, _ -> Some (!index) - | _, _ -> i in - Array.fold_left empty None arr - - (* returns the the location of the nth non-zero - * element in the matrix. Scans column wise. So the nth non-zero element is - * the FIRST non-zero element in the nth non-zero column *) - let nth_nz_location (m: matrix) (n: int): (int*int) option = - let ((n,p),m') = m in - let rec check_col (to_skip: int) (j: int) = - if j <= p then - let (_,col) = get_column m j in - match zero col with - | None -> check_col to_skip (j + 1) - | Some i -> - if to_skip = 0 then - Some (i,j) - else (* we want a later column *) - check_col (to_skip - 1) (j + 1) - else None in - check_col (n - 1) 1 - - (* returns the the location of the first - * non-zero and non-one elt. Scans column wise, from - * left to right. Basically, it ignores columns - * that are all zero or that *) - let fst_nz_no_loc (m: matrix): (int*int) option = - let ((n,p),m') = m in - let rec check_col (j: int) = - if j <= p then - let (_,col) = get_column m j in - match zero col with - | None -> check_col (j + 1) - | Some i -> - match C.compare col.(i-1) C.one with - | Equal -> check_col (j + 1) - | _ -> Some (i,j) - else None in - check_col 1 - - (* Compares two elements in an elt array and returns the greater and its - * index. Is a helper function for find_max_col_index *) - let compare_helper (e1: elt) (e2: elt) (ind1: int) (ind2: int) : (elt*int) = - match C.compare e1 e2 with - | Equal -> (e2, ind2) - | Greater -> (e1, ind1) - | Less -> (e2, ind2) - - (* Finds the element with the greatest absolute value in a column. Is not - * 0-indexed. If two elements are both the maximum value, returns the one with - * the lowest index. Returns None if this element is zero (if column is all 0) - *) - let find_max_col_index (array1: elt array) (start_index: int) : int option = - let rec find_index (max_index: int) (curr_max: elt) (curr_index: int) - (arr: elt array) = - if curr_index = Array.length arr then - (if curr_max = C.zero then None - else Some (max_index+1)) (* Arrays are 0-indexed but matrices aren't *) - else - (match C.compare arr.(curr_index) C.zero with - | Equal -> find_index max_index curr_max (curr_index+1) arr - | Greater -> - (let (el, index) = compare_helper (arr.(curr_index)) - curr_max curr_index max_index in - find_index index el (curr_index+1) arr) - | Less -> - (let abs_curr_elt = C.subtract C.zero arr.(curr_index) in - let (el, index) = compare_helper abs_curr_elt curr_max curr_index - max_index in - find_index index el (curr_index+1) arr)) - in - find_index 0 C.zero (start_index -1) array1 - - (* Basic row operations *) - (* Scales a row by sc *) - let scale_row (m: matrix) (num: int) (sc: elt) : unit = - let (len, row) = get_row m num in - let new_row = Array.map (C.multiply sc) row in - set_row m num new_row - - (* Swaps two rows of a matrix *) - let swap_row (m: matrix) (r1: int) (r2: int) : unit = - let (len1, row1) = get_row m r1 in - let (len2, row2) = get_row m r2 in - let _ = assert (len1 = len2) in - let _ = set_row m r1 row2 in - let _ = set_row m r2 row1 in - () - - (* Subtracts a multiple of r2 from r1 *) - let sub_mult (m: matrix) (r1: int) (r2: int) (sc: elt) : unit = - let (len1, row1) = get_row m r1 in - let (len2, row2) = get_row m r2 in - let _ = assert (len1 = len2) in - for i = 0 to len1 - 1 do (* Arrays are 0-indexed *) - row1.(i) <- C.subtract row1.(i) (C.multiply sc row2.(i)) - done; - set_row m r1 row1 - - (*************** End Helper Functions for Row Reduce ***************) - - (* Returns the row reduced form of a matrix. Is not done in place, but creates - * a new matrix *) - let row_reduce (mat: matrix) : matrix = - let rec row_reduce_h (n_row: int) (n_col: int) (mat2: matrix) : unit = - let ((num_row, num_col), arr) = mat2 in - if (n_col = num_row + 1) then () - else - let (_,col) = get_column mat2 n_col in - match find_max_col_index col n_row with - | None (* Column all 0s *) -> row_reduce_h n_row (n_col+1) mat2 - | Some index -> - begin - swap_row mat2 index n_row; - let pivot = get_elt mat2 (n_row, n_col) in - scale_row mat2 (n_row) (C.divide C.one pivot); - for i = 1 to num_row do - if i <> n_row then sub_mult mat2 i n_row (get_elt mat2 (i,n_col)) - done; - row_reduce_h (n_row+1) (n_col+1) mat2 - end - in - (* Copies the matrix *) - let ((n,p),m) = mat in - let (dim,mat_cp) = empty n p in - for i = 0 to n - 1 do - for j = 0 to p - 1 do - mat_cp.(i).(j) <- m.(i).(j) - done; - done; - let _ = row_reduce_h 1 1 (dim,mat_cp) in (dim,mat_cp) - - (* Pretty prints a matrix: mostly used for testing *) - let print (m: matrix) : unit = - let ((row,col), m') = m in - let pretty_print (_: int) (j: int) (e: elt) = - if j = 1 then - print_string "|" - else (); - C.print e; - print_char ' '; - if j = col then - print_string "|\n" - else () in - iteri pretty_print m - - (*************** End Main Functions ***************) - - (************** Input and Output Functions **********) - let rec read_data (chan: in_channel) : elt list list = - try - let row = input_line chan in - let chars = Helpers.explode row "," in - (List.map C.from_string chars)::read_data chan - with End_of_file -> [] - - (* Load a matrix from a file with filename s *) - let load (s: string): matrix = - try - let chan = open_in s in - let m_list = read_data chan in - let matrix = from_list m_list in - close_in chan; matrix - with - | e -> raise e - - let write_data (((row,col),m): matrix) : string = - let buffer = ref "" in - let to_string (_: int) (j: int) (e: elt) = - buffer := !buffer ^ C.to_string e; - if j = col then - buffer := !buffer ^ "\n" - else - buffer := !buffer ^ "," - in - iteri to_string ((row,col),m); - !buffer - - (* Dumps the matrix to a file with filename name *) - let dump (name: string) (m: matrix) : unit = - try - let outchan = open_out name in - let s = write_data m in - output_string outchan s; - close_out outchan - with - | Sys_error e -> print_string e - - - - - (*************** Optional module functions ***************) - - (* calculates the trace of a matrix *) - let trace (((n,p),m): matrix) : elt = - let rec build (elt: elt) (i: int) = - if i > -1 then - build (C.add m.(i).(i) elt) (i - 1) - else - elt in - if n = p then build C.zero (n - 1) - else raise ImproperDimensions - - (* calculates the transpose of a matrix and retuns a new one *) - let transpose (((n,p),m): matrix) = - let (dim,m') = empty p n in - for i = 0 to n - 1 do - for j = 0 to p - 1 do - m'.(j).(i) <- m.(i).(j) - done; - done; - assert(dim = (p,n)); - ((p,n),m') - - (* Returns the inverse of a matrix. Uses a pretty simple algorithm *) - let inverse (mat: matrix) : matrix = - let ((n,p),m) = mat in - if n = p then - (* create augmented matrix *) - let augmented = empty n (2*n) in - for i = 1 to n do - let (dim,col) = get_column mat i in - let arr = Array.make n C.zero in - begin - assert(dim = n); - arr.(i-1) <- C.one; - set_column augmented i col; - set_column augmented (n + i) arr - end - done; - let augmented' = row_reduce augmented in - (* create the inverted matrix and fill in with appropriate values *) - let inverse = empty n n in - for i = 1 to n do - let (dim, col) = get_column augmented' (n + i) in - let _ = assert(dim = n) in - let _ = set_column inverse i col in - () - done; - inverse - else - raise NonSquare - - (* following our invarient, converts a string into a matrix *) - let from_string (s: string) : matrix = - let convert (row: string) : elt list = - let chars = Helpers.explode row "," in - List.map C.from_string chars in - let rows = Helpers.explode s "|" in - let char_list = List.map convert rows in - from_list char_list - - (***************** HELPER FUNCTIONS FOR DETERMINANT *****************) - (* creates an identity matrix of size n*) - let create_identity (n:int) : matrix = - let (dim,m) = empty n n in - for i = 0 to n - 1 do - m.(i).(i) <- C.one - done; - (dim,m) - - (* Finds the index of the maximum value of an array *) - let find_max_index (arr: elt array) (start_index : int) : int = - let rec find_index (max_index: int) (curr_index: int) = - if curr_index = Array.length arr then max_index+1 - else - match C.compare arr.(curr_index) arr.(max_index) with - | Equal | Less -> find_index max_index (curr_index + 1) - | Greater -> find_index curr_index (curr_index + 1) in - find_index (start_index - 1) start_index - - (* Creates the pivoting matrix for A. Returns swqps. Adapted from - * http://rosettacode.org/wiki/LU_decomposition#Common_Lisp *) - let pivotize (((n,p),m): matrix) : matrix * int = - if n = p then - let swaps = ref 0 in - let pivot_mat = create_identity n in - for j = 1 to n do - let (_,col) = get_column ((n,p),m) j in - let max_index = find_max_index col j in - if max_index <> j then - (swaps := !swaps + 1; swap_row pivot_mat max_index j) - else () - done; - (pivot_mat,!swaps) - else raise ImproperDimensions - - (* decomposes a matrix into a lower triangualar, upper triangualar - * and a pivot matrix. It returns (L,U,P). Adapted from - * http://rosettacode.org/wiki/LU_decomposition#Common_Lisp *) - let lu_decomposition (((n,p),m): matrix) : (matrix*matrix*matrix)*int = - if n = p then - let mat = ((n,p),m) in - let lower, upper, (pivot,s) = empty n n, empty n n, pivotize mat in - let ((x1,y1),l),((x2,y2),u),((x3,y3),p) = lower,upper,pivot in - let ((x,y),mat') = mult pivot mat in - for j = 0 to n - 1 do - l.(j).(j) <- C.one; - for i = 0 to j do - let sum = ref C.zero in - for k = 0 to i - 1 do - sum := C.add (!sum) (C.multiply u.(k).(j) l.(i).(k)) - done; - u.(i).(j) <- C.subtract mat'.(i).(j) (!sum) - done; - for i = j to n - 1 do - let sum = ref C.zero in - for k = 0 to j - 1 do - sum := C.add (!sum) (C.multiply u.(k).(j) l.(i).(k)) - done; - let sub = C.subtract mat'.(i).(j) (!sum) in - l.(i).(j) <- C.divide sub u.(j).(j) - done; - done; - (lower,upper,pivot),s - else raise ImproperDimensions - - (* Computes the determinant of a matrix *) - let determinant (m: matrix) : elt = - try - let ((n,p),m') = m in - if n = p then - let rec triangualar_det (a,mat) curr_index acc = - if curr_index < n then - let acc' = C.multiply mat.(curr_index).(curr_index) acc in - triangualar_det (a,mat) (curr_index + 1) acc' - else acc in - let ((dim1,l),(dim2,u),(dim3,p)),s = lu_decomposition m in - let det1, det2 = triangualar_det (dim1,l) 0 C.one, - triangualar_det (dim2,u) 0 C.one in - if s mod 2 = 0 then C.multiply det1 det2 - else C.subtract C.zero (C.multiply det1 det2) - else raise ImproperDimensions - with - | Failure "create_ratio infinite or undefined rational number" -> C.zero - - - (*************** Optional module functions ***************) - - (*************** Tests ***************) - let print_loc i j = print_string ("Failed @ " ^ string_of_int i ^ " and " ^ - string_of_int j) - - (* Generates a linearly independent matrix *) - let gen_non_independent x y = - let mat = empty x y in - for i = 1 to x-1 do - for j = 1 to y do - set_elt mat (i,j) (C.generate_random (float (x + y)) ()) - done; - done; - let (_,row) = get_row mat 1 in - let _ = set_row mat x row in - mat - - let rec test_empty (times: int) : unit = - if times = 0 then () - else - let dimx, dimy = Random.int times + 1, Random.int times + 1 in - let ((x,y),t_mat) = empty dimx dimy in - assert(x = dimx); - assert(y = dimy); - for i = 0 to dimx - 1 do - for j = 0 to dimy - 1 do - match C.compare C.zero t_mat.(i).(j) with - | Equal -> () - | _ -> print_loc i j - done; - done; - test_empty (times - 1) - - let rec test_map (times: int) : unit = - let rec test_from_0 (index: int) (n: elt) (((x,y),m): matrix) : unit = - if times = index then () - else - let ((x',y'),m') = map (C.add n) ((x,y),m) in - assert(x' = x); - assert(y' = y); - for i = 1 to x do - for j = 1 to y do - let elt = C.add n m.(i-1).(j-1) in - match C.compare elt m'.(i-1).(j-1) with - | Equal -> () - | _ -> print_loc i j - done; - done; - test_from_0 (index + 1) (C.add n C.one) ((x',y'),m') in - let dimx, dimy = Random.int times + 1, Random.int times + 1 in - let m = empty dimx dimy in - test_from_0 0 C.zero m - - - let rec test_get_row (times: int): unit = - if times = 0 then () - else - let dimx, dimy = Random.int times + 1, Random.int times + 1 in - let (extra,t_mat) = map (fun _ -> C.generate_random (float times) ()) - (empty dimx dimy) in - for i = 1 to dimx do - let (dim,row) = get_row (extra, t_mat) i in - assert(dim = dimy); - for j = 1 to dimy do - match C.compare t_mat.(i-1).(j-1) row.(j-1) with - | Equal -> () - | _ -> print_loc i j - done; - done; - test_get_row (times - 1) - - let rec test_get_column (times: int): unit = - if times = 0 then () - else - let dimx, dimy = Random.int times + 1, Random.int times + 1 in - let (extra,t_mat) = map (fun _ -> C.generate_random (float times) ()) - (empty dimx dimy) in - for j = 1 to dimy do - let (dim,column) = get_column (extra, t_mat) j in - assert(dim = dimx); - for i = 1 to dimx do - match C.compare t_mat.(i-1).(j-1) column.(i-1) with - | Equal -> () - | _ -> print_loc i j - done; - done; - test_get_column (times - 1) - - let rec test_get_elt (times: int): unit = - if times = 0 then () - else - let dimx, dimy = Random.int times + 1, Random.int times + 1 in - let (extra,t_mat) = map (fun _ -> C.generate_random (float times) ()) - (empty dimx dimy) in - for i = 1 to dimx do - for j = 1 to dimy do - let elt = get_elt (extra, t_mat) (i,j) in - match C.compare elt t_mat.(i-1).(j-1) with - | Equal -> () - | _ -> print_loc i j - done; - done; - test_get_row (times - 1) - - let rec test_set_row (times: int): unit = - if times = 0 then () - else - let dimx, dimy = Random.int times + 1, Random.int times + 1 in - let (extra,t_mat) = map (fun _ -> C.generate_random (float times) ()) - (empty dimx dimy) in - let rand_array = Array.map (fun _ -> C.generate_random (float times) ()) - (Array.make dimy C.one) in - for i = 1 to dimx do - let _ = set_row (extra,t_mat) i rand_array in - for j = 1 to dimy do - match C.compare t_mat.(i-1).(j-1) rand_array.(j-1) with - | Equal -> () - | _ -> print_loc i j - done; - done; - test_set_row (times - 1) - - - let rec test_set_column (times: int): unit = - if times = 0 then () - else - let dimx, dimy = Random.int times + 1, Random.int times + 1 in - let (extra,t_mat) = map (fun _ -> C.generate_random (float times) ()) - (empty dimx dimy) in - let rand_array = Array.map (fun _ -> C.generate_random (float times) ()) - (Array.make dimx C.one) in - for j = 1 to dimy do - let _ = set_column (extra,t_mat) j rand_array in - for i = 1 to dimx do - match C.compare t_mat.(i-1).(j-1) rand_array.(i-1) with - | Equal -> () - | _ -> print_loc i j - done; - done; - test_set_row (times - 1) - - let rec test_reduce (times: int) : unit = - if times = 0 then () - else - let dimx, dimy = Random.int times + 1, Random.int times + 1 in - let t_mat = map (C.add C.one) (empty dimx dimy) in - let result = reduce C.add C.zero t_mat in - let e_result = C.generate_x (float (dimx * dimy)) () in - assert(e_result = result); - test_reduce (times - 1) - - let rec test_mult (times: int) : unit = - if times = 0 then () - else - let dimx = Random.int times + 1 in - let identity = create_identity dimx in - let t_mat = map (fun x -> C.generate_random (float times) ()) - (empty dimx dimx) in - let result = mult identity t_mat in - for i = 1 to dimx do - for j = 1 to dimx do - assert(get_elt t_mat (i,j) = get_elt result (i,j)) - done; - done; - test_mult (times - 1) - - let rec test_add (times: int) : unit = - if times = 0 then () - else - let dimx= Random.int times + 1 in - let identity = create_identity dimx in - let t_mat = map (fun x -> C.generate_random (float times) ()) - (empty dimx dimx) in - let result = add identity t_mat in - for i = 1 to dimx do - for j = 1 to dimx do - if i <> j then assert(get_elt t_mat (i,j) = get_elt result (i,j)) - else assert((C.add C.one (get_elt t_mat (i,j))) = get_elt result (i,j)) - done; - done; - test_add (times - 1) - - let rec test_scale (times: int) : unit = - if times = 0 then () - else - let dimx,dimy = Random.int times + 1, Random.int times + 1 in - let scalar = C.generate_random (float times) () in - let t_mat = map (fun x -> C.generate_random (float times) ()) - (empty dimx dimy) in - let result = scale t_mat scalar in - for i = 1 to dimx do - for j = 1 to dimy do - assert((C.multiply scalar (get_elt t_mat (i,j))) = get_elt result (i,j)) - done; - done; - test_add (times - 1) - - let rec test_row_reduce (times: int) : unit = - let dimx = Random.int times + 2 in - let independent = gen_non_independent dimx dimx in - let identity = create_identity dimx in - let reduced = row_reduce independent in - try - assert (not(reduced = identity)) - with Assert_failure _ -> (print independent;print reduced;print identity) - - let rec test_determinamnt (times: int) : unit = - let dimx = Random.int times + 2 in - let independent = gen_non_independent dimx dimx in - let det = determinant independent in - assert(det = C.zero) - - - - (*************** End Test Functions ***************) - - let run_tests times = - test_empty times; - test_map times; - test_get_row times; - test_get_column times; - test_get_elt times; - test_set_row times; - test_set_column times; - test_reduce times; - test_add times; - test_mult times; - test_scale times; - test_row_reduce times; - test_determinamnt times; - () - -end - -(* Creating the Matrix Library module! *) -module EltMatrix = MakeMatrix(Elts) +open Order +open Elts + +(* Functor so we can Abstract away! *) +module MakeMatrix (C: EltsI.ORDERED_AND_OPERATIONAL) : + (MatrixI.MATRIX with type elt = C.t) = +struct + + (*************** Exceptions ***************) + + exception NonSquare + exception ImproperDimensions + + (*************** End Exceptions ***************) + + (*************** Types ***************) + + type elt = C.t + + (* A matrix is a pair of dimension (n x p) and a array of arrays + * the first array is the row (n) and the second the column (p) *) + type matrix = (int * int) * (elt array array) + + (*************** End Types ***************) + + (*************** Base Functions ***************) + + (* catching negative dimensions AND 0 dimensions and too large + * of a dimension so we don't have to worry about it later *) + let empty (rows: int) (columns: int) : matrix = + if rows > 0 && columns > 0 then + try + let m = Array.make_matrix rows columns C.zero in ((rows,columns),m) + with Invalid_argument "index out of bounds" -> + raise ImproperDimensions + else (* dimension is negative or 0 *) + raise ImproperDimensions + + (*************** End Base Functions ***************) + + (*************** Helper Functions ***************) + + (* get's the nth row of a matrix and returns (r, row) where r is the length + * of the row and row is a COPY of the original row. For example, calling + * calling get_row m 1 will return (3, |1 3 4 |) + * ________ + * m = | 1 3 4 | + * |*2 5 6 | + *) + (* aside: we don't check whether n < 1 because of our matrix invariant *) + let get_row (((n,p),m): matrix) (row: int) : int * elt array = + if row <= n then + let row' = Array.map (fun x -> x) m.(row - 1) in + (p, row') + else + raise (Failure "Row out of bounds.") + + (* similar to get_row. For m, get_column m 1 will return (2, |1 2|) *) + let get_column (((n,p),m): matrix) (column: int) : int * elt array = + if column <= p then + begin + let column' = Array.make n C.zero in + for i = 0 to n - 1 do + column'.(i) <- m.(i).(column - 1) + done; + (n, column') + end + else + raise (Failure "Column out of bounds.") + + (* sets the nth row of the matrix m to the specified array a. + * This is done IN-PLACE. Therefore the function returns unit. You should + * nonetheless enfore immutability whenever possible. For a clarification on + * what nth row means, look at comment for get_row above. *) + let set_row (((n,p),m): matrix) (row: int) (a: elt array) : unit = + if row <= n then + begin + assert(Array.length a = p); + for i = 0 to p - 1 do + m.(row - 1).(i) <- a.(i) + done; + end + else + raise (Failure "Row out of bounds.") + + (* Similar to set_row but sets the nth column instead *) + let set_column (((n,p),m): matrix) (column: int) (a: elt array) : unit = + if column <= p then + begin + assert(Array.length a = n); + for i = 0 to n - 1 do + m.(i).(column - 1) <- a.(i) + done; + end + else + raise (Failure "Column out of bounds.") + + (* returns the ij-th element of a matrix (not-zero indexed) *) + let get_elt (((n,p),m): matrix) ((i,j): int*int) : elt = + if i <= n && j <= p then + m.(i - 1).(j - 1) + else + raise ImproperDimensions + + (* Changes the i,j-th element of a matrix to e. Is not zero-indexed, and + * changes the matrix in place *) + let set_elt (((n,p),m): matrix) ((i,j): int*int) (e: elt) : unit = + if i <= n && j <= p then + m.(i - 1).(j - 1) <- e + else + raise ImproperDimensions + + (* similar to map, but applies to function to the entire matrix + * Returns a new matrix *) + let map (f: elt -> elt) (mat: matrix) : matrix = + let (dim,m) = mat in + (dim, Array.map (Array.map f) m) + + (* Just some wrapping of Array.iter made for Matrices! *) + let iter (f: elt -> unit) (mat: matrix) : unit = + let dim,m = mat in + Array.iter (Array.iter f) m + + (* Just some wrapping of Array.iteri. Useful for pretty + * printing matrix. The index is (i,j). NOT zero-indexed *) + let iteri (f: int -> int -> elt -> unit) (mat: matrix) : unit = + let dim,m = mat in + Array.iteri (fun i row -> Array.iteri (fun j e -> f (i+1) (j+1) e) row) m + + (* folds over each row using base case u and function f *) + (* could be a bit more efficient? *) + let reduce (f: 'a -> elt -> 'a) (u: 'a) (((p,q),m): matrix) : 'a = + let total = ref u in + for i = 0 to p - 1 do + for j = 0 to q - 1 do + total := f (!total) m.(i).(j) + done; + done; + !total + + (* given two arrays, this will calculate their dot product *) + (* It seems a little funky, but this is done for efficiency's sake. + * In short, it tries to multiply each element by it's respective + * element until the one array is indexed out of bounds. If the + * other array is also out of bounds, then it returns their value. + * Otherwise, the arrays were the wrong size and raises ImproperDimension + + THE ABOVE COMMENT HAS NOT BEEN IMPLEMENTED + + Instead we calculate the length before starting + *) + let dot (v1: elt array) (v2: elt array) : elt = + let rec dotting (i: int) (total: elt) : elt = + if i = 0 then total + else + let curr = C.multiply v1.(i-1) v2.(i-1) in + dotting (i - 1) (C.add curr total) in + let len1, len2 = Array.length v1, Array.length v2 in + if len1 = len2 then dotting len1 C.zero + else raise ImproperDimensions + + (* function to expose the dimensions of a matrix *) + let get_dimensions (m: matrix) : (int * int) = + let ((x,y),m') = m in (x,y) + + (*************** End Helper Functions ***************) + + + (*************** Primary Matrix Functions ***************) + + (* scales a matrix by the appropriate factor *) + let scale (m: matrix) (sc: elt) : matrix = map (C.multiply sc) m + + (* Generates a matrix from a list of lists. The inners lists are the rows *) + let from_list (lsts : elt list list) : matrix = + let rec check_length (length: int) (lst: elt list) : int = + if List.length lst = length then length + else raise ImproperDimensions in + let p = List.length lsts in + match lsts with + | [] -> raise ImproperDimensions + | hd::tl -> + let len = List.length hd in + if List.fold_left check_length len tl = len then + ((p,len),Array.map Array.of_list (Array.of_list lsts)) + else + raise ImproperDimensions + + (* Adds two matrices. They must have the same dimensions *) + let add ((dim1,m1): matrix) ((dim2,m2): matrix) : matrix = + if dim1 = dim2 then + let n, p = dim1 in + let (dim', sum_m) = empty n p in + for i = 0 to n - 1 do + for j = 0 to p - 1 do + sum_m.(i).(j) <- C.add m1.(i).(j) m2.(i).(j) + done; + done; + (dim',sum_m) + else + raise ImproperDimensions + + + (* Multiplies two matrices. If the matrices have dimensions m x n and p x q, n + * and p must be equal, and the resulting matrix will have dimension n x q *) + let mult (matrix1: matrix) (matrix2: matrix) : matrix = + let ((m,n),m1), ((p,q),m2) = matrix1, matrix2 in + if n = p then + let (dim, result) = empty m q in + for i = 0 to m - 1 do + for j = 0 to q - 1 do + let (_,row), (_,column) = get_row matrix1 (i + 1), + get_column matrix2 (j + 1) in + result.(i).(j) <- dot row column + done; + done; + (dim,result) + else + raise ImproperDimensions + + (*************** Helper Functions for Row Reduce ***************) + + (* returns the index of the first non-zero elt in an array*) + let zero (arr: elt array) : int option = + let index = ref 1 in + let empty (i: int option) (e: elt) : int option = + match i, C.compare e C.zero with + | None, Equal -> (index := !index + 1; None) + | None, _ -> Some (!index) + | _, _ -> i in + Array.fold_left empty None arr + + (* returns the the location of the nth non-zero + * element in the matrix. Scans column wise. So the nth non-zero element is + * the FIRST non-zero element in the nth non-zero column *) + let nth_nz_location (m: matrix) (n: int): (int*int) option = + let ((n,p),m') = m in + let rec check_col (to_skip: int) (j: int) = + if j <= p then + let (_,col) = get_column m j in + match zero col with + | None -> check_col to_skip (j + 1) + | Some i -> + if to_skip = 0 then + Some (i,j) + else (* we want a later column *) + check_col (to_skip - 1) (j + 1) + else None in + check_col (n - 1) 1 + + (* returns the the location of the first + * non-zero and non-one elt. Scans column wise, from + * left to right. Basically, it ignores columns + * that are all zero or that *) + let fst_nz_no_loc (m: matrix): (int*int) option = + let ((n,p),m') = m in + let rec check_col (j: int) = + if j <= p then + let (_,col) = get_column m j in + match zero col with + | None -> check_col (j + 1) + | Some i -> + match C.compare col.(i-1) C.one with + | Equal -> check_col (j + 1) + | _ -> Some (i,j) + else None in + check_col 1 + + (* Compares two elements in an elt array and returns the greater and its + * index. Is a helper function for find_max_col_index *) + let compare_helper (e1: elt) (e2: elt) (ind1: int) (ind2: int) : (elt*int) = + match C.compare e1 e2 with + | Equal -> (e2, ind2) + | Greater -> (e1, ind1) + | Less -> (e2, ind2) + + (* Finds the element with the greatest absolute value in a column. Is not + * 0-indexed. If two elements are both the maximum value, returns the one with + * the lowest index. Returns None if this element is zero (if column is all 0) + *) + let find_max_col_index (array1: elt array) (start_index: int) : int option = + let rec find_index (max_index: int) (curr_max: elt) (curr_index: int) + (arr: elt array) = + if curr_index = Array.length arr then + (if curr_max = C.zero then None + else Some (max_index+1)) (* Arrays are 0-indexed but matrices aren't *) + else + (match C.compare arr.(curr_index) C.zero with + | Equal -> find_index max_index curr_max (curr_index+1) arr + | Greater -> + (let (el, index) = compare_helper (arr.(curr_index)) + curr_max curr_index max_index in + find_index index el (curr_index+1) arr) + | Less -> + (let abs_curr_elt = C.subtract C.zero arr.(curr_index) in + let (el, index) = compare_helper abs_curr_elt curr_max curr_index + max_index in + find_index index el (curr_index+1) arr)) + in + find_index 0 C.zero (start_index -1) array1 + + (* Basic row operations *) + (* Scales a row by sc *) + let scale_row (m: matrix) (num: int) (sc: elt) : unit = + let (len, row) = get_row m num in + let new_row = Array.map (C.multiply sc) row in + set_row m num new_row + + (* Swaps two rows of a matrix *) + let swap_row (m: matrix) (r1: int) (r2: int) : unit = + let (len1, row1) = get_row m r1 in + let (len2, row2) = get_row m r2 in + let _ = assert (len1 = len2) in + let _ = set_row m r1 row2 in + let _ = set_row m r2 row1 in + () + + (* Subtracts a multiple of r2 from r1 *) + let sub_mult (m: matrix) (r1: int) (r2: int) (sc: elt) : unit = + let (len1, row1) = get_row m r1 in + let (len2, row2) = get_row m r2 in + let _ = assert (len1 = len2) in + for i = 0 to len1 - 1 do (* Arrays are 0-indexed *) + row1.(i) <- C.subtract row1.(i) (C.multiply sc row2.(i)) + done; + set_row m r1 row1 + + (*************** End Helper Functions for Row Reduce ***************) + + (* Returns the row reduced form of a matrix. Is not done in place, but creates + * a new matrix *) + let row_reduce (mat: matrix) : matrix = + let rec row_reduce_h (n_row: int) (n_col: int) (mat2: matrix) : unit = + let ((num_row, num_col), arr) = mat2 in + if (n_col = num_row + 1) then () + else + let (_,col) = get_column mat2 n_col in + match find_max_col_index col n_row with + | None (* Column all 0s *) -> row_reduce_h n_row (n_col+1) mat2 + | Some index -> + begin + swap_row mat2 index n_row; + let pivot = get_elt mat2 (n_row, n_col) in + scale_row mat2 (n_row) (C.divide C.one pivot); + for i = 1 to num_row do + if i <> n_row then sub_mult mat2 i n_row (get_elt mat2 (i,n_col)) + done; + row_reduce_h (n_row+1) (n_col+1) mat2 + end + in + (* Copies the matrix *) + let ((n,p),m) = mat in + let (dim,mat_cp) = empty n p in + for i = 0 to n - 1 do + for j = 0 to p - 1 do + mat_cp.(i).(j) <- m.(i).(j) + done; + done; + let _ = row_reduce_h 1 1 (dim,mat_cp) in (dim,mat_cp) + + (* Pretty prints a matrix: mostly used for testing *) + let print (m: matrix) : unit = + let ((row,col), m') = m in + let pretty_print (_: int) (j: int) (e: elt) = + if j = 1 then + print_string "|" + else (); + C.print e; + print_char ' '; + if j = col then + print_string "|\n" + else () in + iteri pretty_print m + + (*************** End Main Functions ***************) + + (************** Input and Output Functions **********) + let rec read_data (chan: in_channel) : elt list list = + try + let row = input_line chan in + let chars = Helpers.explode row "," in + (List.map C.from_string chars)::read_data chan + with End_of_file -> [] + + (* Load a matrix from a file with filename s *) + let load (s: string): matrix = + try + let chan = open_in s in + let m_list = read_data chan in + let matrix = from_list m_list in + close_in chan; matrix + with + | e -> raise e + + let write_data (((row,col),m): matrix) : string = + let buffer = ref "" in + let to_string (_: int) (j: int) (e: elt) = + buffer := !buffer ^ C.to_string e; + if j = col then + buffer := !buffer ^ "\n" + else + buffer := !buffer ^ "," + in + iteri to_string ((row,col),m); + !buffer + + (* Dumps the matrix to a file with filename name *) + let dump (name: string) (m: matrix) : unit = + try + let outchan = open_out name in + let s = write_data m in + output_string outchan s; + close_out outchan + with + | Sys_error e -> print_string e + + + + + (*************** Optional module functions ***************) + + (* calculates the trace of a matrix *) + let trace (((n,p),m): matrix) : elt = + let rec build (elt: elt) (i: int) = + if i > -1 then + build (C.add m.(i).(i) elt) (i - 1) + else + elt in + if n = p then build C.zero (n - 1) + else raise ImproperDimensions + + (* calculates the transpose of a matrix and retuns a new one *) + let transpose (((n,p),m): matrix) = + let (dim,m') = empty p n in + for i = 0 to n - 1 do + for j = 0 to p - 1 do + m'.(j).(i) <- m.(i).(j) + done; + done; + assert(dim = (p,n)); + ((p,n),m') + + (* Returns the inverse of a matrix. Uses a pretty simple algorithm *) + let inverse (mat: matrix) : matrix = + let ((n,p),m) = mat in + if n = p then + (* create augmented matrix *) + let augmented = empty n (2*n) in + for i = 1 to n do + let (dim,col) = get_column mat i in + let arr = Array.make n C.zero in + begin + assert(dim = n); + arr.(i-1) <- C.one; + set_column augmented i col; + set_column augmented (n + i) arr + end + done; + let augmented' = row_reduce augmented in + (* create the inverted matrix and fill in with appropriate values *) + let inverse = empty n n in + for i = 1 to n do + let (dim, col) = get_column augmented' (n + i) in + let _ = assert(dim = n) in + let _ = set_column inverse i col in + () + done; + inverse + else + raise NonSquare + + (* following our invarient, converts a string into a matrix *) + let from_string (s: string) : matrix = + let convert (row: string) : elt list = + let chars = Helpers.explode row "," in + List.map C.from_string chars in + let rows = Helpers.explode s "|" in + let char_list = List.map convert rows in + from_list char_list + + (***************** HELPER FUNCTIONS FOR DETERMINANT *****************) + (* creates an identity matrix of size n*) + let create_identity (n:int) : matrix = + let (dim,m) = empty n n in + for i = 0 to n - 1 do + m.(i).(i) <- C.one + done; + (dim,m) + + (* Finds the index of the maximum value of an array *) + let find_max_index (arr: elt array) (start_index : int) : int = + let rec find_index (max_index: int) (curr_index: int) = + if curr_index = Array.length arr then max_index+1 + else + match C.compare arr.(curr_index) arr.(max_index) with + | Equal | Less -> find_index max_index (curr_index + 1) + | Greater -> find_index curr_index (curr_index + 1) in + find_index (start_index - 1) start_index + + (* Creates the pivoting matrix for A. Returns swqps. Adapted from + * http://rosettacode.org/wiki/LU_decomposition#Common_Lisp *) + let pivotize (((n,p),m): matrix) : matrix * int = + if n = p then + let swaps = ref 0 in + let pivot_mat = create_identity n in + for j = 1 to n do + let (_,col) = get_column ((n,p),m) j in + let max_index = find_max_index col j in + if max_index <> j then + (swaps := !swaps + 1; swap_row pivot_mat max_index j) + else () + done; + (pivot_mat,!swaps) + else raise ImproperDimensions + + (* decomposes a matrix into a lower triangualar, upper triangualar + * and a pivot matrix. It returns (L,U,P). Adapted from + * http://rosettacode.org/wiki/LU_decomposition#Common_Lisp *) + let lu_decomposition (((n,p),m): matrix) : (matrix*matrix*matrix)*int = + if n = p then + let mat = ((n,p),m) in + let lower, upper, (pivot,s) = empty n n, empty n n, pivotize mat in + let ((x1,y1),l),((x2,y2),u),((x3,y3),p) = lower,upper,pivot in + let ((x,y),mat') = mult pivot mat in + for j = 0 to n - 1 do + l.(j).(j) <- C.one; + for i = 0 to j do + let sum = ref C.zero in + for k = 0 to i - 1 do + sum := C.add (!sum) (C.multiply u.(k).(j) l.(i).(k)) + done; + u.(i).(j) <- C.subtract mat'.(i).(j) (!sum) + done; + for i = j to n - 1 do + let sum = ref C.zero in + for k = 0 to j - 1 do + sum := C.add (!sum) (C.multiply u.(k).(j) l.(i).(k)) + done; + let sub = C.subtract mat'.(i).(j) (!sum) in + l.(i).(j) <- C.divide sub u.(j).(j) + done; + done; + (lower,upper,pivot),s + else raise ImproperDimensions + + (* Computes the determinant of a matrix *) + let determinant (m: matrix) : elt = + try + let ((n,p),m') = m in + if n = p then + let rec triangualar_det (a,mat) curr_index acc = + if curr_index < n then + let acc' = C.multiply mat.(curr_index).(curr_index) acc in + triangualar_det (a,mat) (curr_index + 1) acc' + else acc in + let ((dim1,l),(dim2,u),(dim3,p)),s = lu_decomposition m in + let det1, det2 = triangualar_det (dim1,l) 0 C.one, + triangualar_det (dim2,u) 0 C.one in + if s mod 2 = 0 then C.multiply det1 det2 + else C.subtract C.zero (C.multiply det1 det2) + else raise ImproperDimensions + with + | Failure "create_ratio infinite or undefined rational number" -> C.zero + + + (*************** Optional module functions ***************) + + (*************** Tests ***************) + let print_loc i j = print_string ("Failed @ " ^ string_of_int i ^ " and " ^ + string_of_int j) + + (* Generates a linearly independent matrix *) + let gen_non_independent x y = + let mat = empty x y in + for i = 1 to x-1 do + for j = 1 to y do + set_elt mat (i,j) (C.generate_random (float (x + y)) ()) + done; + done; + let (_,row) = get_row mat 1 in + let _ = set_row mat x row in + mat + + let rec test_empty (times: int) : unit = + if times = 0 then () + else + let dimx, dimy = Random.int times + 1, Random.int times + 1 in + let ((x,y),t_mat) = empty dimx dimy in + assert(x = dimx); + assert(y = dimy); + for i = 0 to dimx - 1 do + for j = 0 to dimy - 1 do + match C.compare C.zero t_mat.(i).(j) with + | Equal -> () + | _ -> print_loc i j + done; + done; + test_empty (times - 1) + + let rec test_map (times: int) : unit = + let rec test_from_0 (index: int) (n: elt) (((x,y),m): matrix) : unit = + if times = index then () + else + let ((x',y'),m') = map (C.add n) ((x,y),m) in + assert(x' = x); + assert(y' = y); + for i = 1 to x do + for j = 1 to y do + let elt = C.add n m.(i-1).(j-1) in + match C.compare elt m'.(i-1).(j-1) with + | Equal -> () + | _ -> print_loc i j + done; + done; + test_from_0 (index + 1) (C.add n C.one) ((x',y'),m') in + let dimx, dimy = Random.int times + 1, Random.int times + 1 in + let m = empty dimx dimy in + test_from_0 0 C.zero m + + + let rec test_get_row (times: int): unit = + if times = 0 then () + else + let dimx, dimy = Random.int times + 1, Random.int times + 1 in + let (extra,t_mat) = map (fun _ -> C.generate_random (float times) ()) + (empty dimx dimy) in + for i = 1 to dimx do + let (dim,row) = get_row (extra, t_mat) i in + assert(dim = dimy); + for j = 1 to dimy do + match C.compare t_mat.(i-1).(j-1) row.(j-1) with + | Equal -> () + | _ -> print_loc i j + done; + done; + test_get_row (times - 1) + + let rec test_get_column (times: int): unit = + if times = 0 then () + else + let dimx, dimy = Random.int times + 1, Random.int times + 1 in + let (extra,t_mat) = map (fun _ -> C.generate_random (float times) ()) + (empty dimx dimy) in + for j = 1 to dimy do + let (dim,column) = get_column (extra, t_mat) j in + assert(dim = dimx); + for i = 1 to dimx do + match C.compare t_mat.(i-1).(j-1) column.(i-1) with + | Equal -> () + | _ -> print_loc i j + done; + done; + test_get_column (times - 1) + + let rec test_get_elt (times: int): unit = + if times = 0 then () + else + let dimx, dimy = Random.int times + 1, Random.int times + 1 in + let (extra,t_mat) = map (fun _ -> C.generate_random (float times) ()) + (empty dimx dimy) in + for i = 1 to dimx do + for j = 1 to dimy do + let elt = get_elt (extra, t_mat) (i,j) in + match C.compare elt t_mat.(i-1).(j-1) with + | Equal -> () + | _ -> print_loc i j + done; + done; + test_get_row (times - 1) + + let rec test_set_row (times: int): unit = + if times = 0 then () + else + let dimx, dimy = Random.int times + 1, Random.int times + 1 in + let (extra,t_mat) = map (fun _ -> C.generate_random (float times) ()) + (empty dimx dimy) in + let rand_array = Array.map (fun _ -> C.generate_random (float times) ()) + (Array.make dimy C.one) in + for i = 1 to dimx do + let _ = set_row (extra,t_mat) i rand_array in + for j = 1 to dimy do + match C.compare t_mat.(i-1).(j-1) rand_array.(j-1) with + | Equal -> () + | _ -> print_loc i j + done; + done; + test_set_row (times - 1) + + + let rec test_set_column (times: int): unit = + if times = 0 then () + else + let dimx, dimy = Random.int times + 1, Random.int times + 1 in + let (extra,t_mat) = map (fun _ -> C.generate_random (float times) ()) + (empty dimx dimy) in + let rand_array = Array.map (fun _ -> C.generate_random (float times) ()) + (Array.make dimx C.one) in + for j = 1 to dimy do + let _ = set_column (extra,t_mat) j rand_array in + for i = 1 to dimx do + match C.compare t_mat.(i-1).(j-1) rand_array.(i-1) with + | Equal -> () + | _ -> print_loc i j + done; + done; + test_set_row (times - 1) + + let rec test_reduce (times: int) : unit = + if times = 0 then () + else + let dimx, dimy = Random.int times + 1, Random.int times + 1 in + let t_mat = map (C.add C.one) (empty dimx dimy) in + let result = reduce C.add C.zero t_mat in + let e_result = C.generate_x (float (dimx * dimy)) () in + assert(e_result = result); + test_reduce (times - 1) + + let rec test_mult (times: int) : unit = + if times = 0 then () + else + let dimx = Random.int times + 1 in + let identity = create_identity dimx in + let t_mat = map (fun x -> C.generate_random (float times) ()) + (empty dimx dimx) in + let result = mult identity t_mat in + for i = 1 to dimx do + for j = 1 to dimx do + assert(get_elt t_mat (i,j) = get_elt result (i,j)) + done; + done; + test_mult (times - 1) + + let rec test_add (times: int) : unit = + if times = 0 then () + else + let dimx= Random.int times + 1 in + let identity = create_identity dimx in + let t_mat = map (fun x -> C.generate_random (float times) ()) + (empty dimx dimx) in + let result = add identity t_mat in + for i = 1 to dimx do + for j = 1 to dimx do + if i <> j then assert(get_elt t_mat (i,j) = get_elt result (i,j)) + else assert((C.add C.one (get_elt t_mat (i,j))) = get_elt result (i,j)) + done; + done; + test_add (times - 1) + + let rec test_scale (times: int) : unit = + if times = 0 then () + else + let dimx,dimy = Random.int times + 1, Random.int times + 1 in + let scalar = C.generate_random (float times) () in + let t_mat = map (fun x -> C.generate_random (float times) ()) + (empty dimx dimy) in + let result = scale t_mat scalar in + for i = 1 to dimx do + for j = 1 to dimy do + assert((C.multiply scalar (get_elt t_mat (i,j))) = get_elt result (i,j)) + done; + done; + test_add (times - 1) + + let rec test_row_reduce (times: int) : unit = + let dimx = Random.int times + 2 in + let independent = gen_non_independent dimx dimx in + let identity = create_identity dimx in + let reduced = row_reduce independent in + try + assert (not(reduced = identity)) + with Assert_failure _ -> (print independent;print reduced;print identity) + + let rec test_determinamnt (times: int) : unit = + let dimx = Random.int times + 2 in + let independent = gen_non_independent dimx dimx in + let det = determinant independent in + assert(det = C.zero) + + + + (*************** End Test Functions ***************) + + let run_tests times = + test_empty times; + test_map times; + test_get_row times; + test_get_column times; + test_get_elt times; + test_set_row times; + test_set_column times; + test_reduce times; + test_add times; + test_mult times; + test_scale times; + test_row_reduce times; + test_determinamnt times; + () + +end + +(* Creating the Matrix Library module! *) +module EltMatrix = MakeMatrix(Elts) diff --git a/info.txt b/info.txt index b8a942c..c380484 100644 --- a/info.txt +++ b/info.txt @@ -1,2 +1,2 @@ -title: Matrix Module and Simplex Algorithm +title: Matrix Module and Simplex Algorithm seas-usernames: andyshi, dzhou, luisperez, zihaowang01 \ No newline at end of file diff --git a/maketop b/maketop deleted file mode 100644 index 6f0a7c7..0000000 --- a/maketop +++ /dev/null @@ -1 +0,0 @@ -ocamlmktop -o simplextop nums.cma Order.cmo Elts.cmo Helpers.cmo Matrix.cmo SimplexI.cmo Interface.cmo diff --git a/report.tex b/report.tex index 7673f7f..8d98586 100644 --- a/report.tex +++ b/report.tex @@ -1,328 +1,328 @@ -\documentclass[letterpaper,12pt]{article} -% Change the header if you're going to change this. - -% Possible packages - uncomment to use -%\usepackage{amsmath} % Needed for math stuff -%\usepackage{lastpage} % Finds last -%\usepackage{amssymb} % Needed for some math symbols -%\usepackage{graphicx} % Needed for graphics -\usepackage[usenames,dvipsnames]{xcolor} % Needed for graphics and color -%\usepackage{setspace} % Needed to set single/double spacing -\usepackage{hyperref} - -%Suggested by TeXworks -\usepackage[utf8]{inputenc} % set input encoding (not needed with XeLaTeX) -\usepackage{verbatim} % adds environment for commenting out blocks of text & for better verbatim - -% Sets the page margins to 1 inch each side -\usepackage[margin=1in]{geometry} - -% Uncomment this section if you wish to have a header. -\usepackage{fancyhdr} -\pagestyle{fancy} -\renewcommand{\headrulewidth}{0.5pt} % customise the layout... -\lhead{CS51 Final Proect Report} \chead{} \rhead{May 5, 2013} -\lfoot{} \cfoot{\thepage} \rfoot{} - -\frenchspacing -\setlength\parindent{0pt} -\setlength\parskip{2ex} - -\newcommand{\annot}[1]{\textbf{\color{BrickRed} [#1]}} - -\newcommand{\vect}[1]{\mathbf{#1}} - -\begin{document} -\title{CS 51 Final Project Report} -\author{ -Luis Perez $|$ luisperez@college.harvard.edu \\ -Andy Shi $|$ andyshi@college.harvard.edu \\ -Zihao Wang $|$ zihaowang01@college.harvard.edu \\ -Ding Zhou $|$ dzhou@college.harvard.edu -} -\date{May 5, 2013} -\maketitle - -Link to demo video: \url{} - -\section{Overview} - -In the real world, people love to make money. Even in virtual worlds (bitcoin?), -people love to make money. Basically, many problems in business and even -everyday life (but especially business), boil down to maximizing profit or -minimizing costs. Of course, when you speak of minimization or maximization -there are always constraints. Resources are finite and time is precious. -Therefore, in an abstract sense, the problem boils down to minimizing or -maximizing some quantity under certain constraints. For example, you might want -to maximize the number of textbooks $x,y,z$ you buy under the constraints that you -can't buy more than five from one company, the total cost can't be more -than \$500, and the last two books you buy must cost more than \$75. This problem -is precisely the type of problem that our program can now solve. How did we do -this? - -To start off, we created a framework. We implemented a matrix module with -standard matrix operations such as multiplication, addition, map (where scaling -follows) as well as many extra features such as row reduction, determinant, -trace, and many more. The matrix module basically allows a user to do almost -anything imaginable with a matrix. In addition, we provided read and write -functionality for easy loading of data, and printing functionality for -straightforward viewing of the abstract matrix. However the intellectually -challenging part of our project is the implementation of the simplex algorithm. -We needed to create a matrix module simply because OCaml doesn't provide one -(also because we just wanted to implement a matrix module, of course). We use -many of the functions in the matrix module to create and manipulate simplex -systems within our simplex module. This allowed us to abstract away yet another -part of our project. - -The simplex algorithm is a type of linear programming that allows one to -optimize a linear equation given linear constraints (like the example problem at -the beginning). In mathematical terms, the simplex program operates on linear -programs in standard form. That is, it seeks to minimize (or maximize) -$C\vect{x}$ subject to $A\vect{x} = \vect{b} $ (Wikipedia). $\vect{x}$ and -$\vect{b}$ are vectors. $x$ contains the variables in the equation we want to -optimize (the objective), and $b$ contains the constants in the constraints. $C$ -and $A$ are matrices, and they represent the coefficients in the objective and -the constraints, respectively. Let's consider the inequalities: $30x + 15y + 18z -\geq 500, x \geq 5,y \geq 5, z \geq 5$, and $15y+18z \geq 75$ and the objective -functions $x +y + z$ where $x, y, z \geq 0$. As you will notice, this is -precisely the system of equations necessary to solve our book example above. The -algorithm converts the inequalities into equalities by adding positive slack -variables. Then it represents this set of equations in the matrix form where the -coefficients of the terms we maximize are on the first row and each column -represents a variable. This provides affordable storage of most of the -information necessary to solve the linear program. Furthermore, we then define -basic and non-basic variables (kept track of in a list), which together with the -aforementioned matrix, create a simplex system. - -The first function necessary for simplex is called \verb@simple_solve@. This -function, in abstract terms, finds the minimum value of a passed in system (as -described above). When we call the function \verb@simple_solve@, the program -finds the entering and leaving variables which are the variable that would -exchange position from basic and non-basic (slack to non-slack in the initial -case). Here, the algorithm is making the assumption that the system has a -feasible solution where all of our non-basic variables are zero, and therefore -exchanges the role of one of these variables with the corresponding basic -variable. In geometric terms, we are moving along one edge of the n-gon formed -by the constraints and attempting to minimize. This step is accomplished by the -helper function pivot. The function \verb@simple_solve@ continues to loop -(moving along edges and minimizing) until the first row---which represents the -term we want to optimize consists of only non-positive values or there is no -possible entering variable (ie, any change in the variables will only increase -our answer, not minimize it). In the first case, algebraically, it means that -the constraint can only increase if the variables change. Therefore, we have -found our minimum value which is the constant stored at the end of the first -row. This is the optimized value of the objective function. In the second case, -where no suitable entering variable can be found, it means that the system can -be decreased infinitely as pleased and therefore it is Unbounded. In the case -where a solution exists, one possible value for the original variables can also -be found by looking at their corresponding columns and checking to see if they -are basic. If there are an infinite number of these, our algorithm returns only -one. - -What if we had the constraint $x + y + z \geq -10$? As stated above, our -\verb@simple_solve@ makes the assumption that one of the feasible solutions is -the zero solution, but if we initially set $x, y$, and $z$ to 0 (our algorithm -generates a ``basic solution'' by doing this step) the constraint does not hold. -We therefore have a helper function called \verb@initialize_simplex@ that -converts everything into the ``simple'' form so we can just recursively call -\verb@simple_solve@. It does this by following a very complex algorithm better -detailed below. Our final simplex algorithm now accounts for all cases as long -as the equations are linear. It also handles not only less than or equal -constraints, but also greater than or equal to and equality constraints. -Finally, not only can it minimize, but it can also maximize. You can follow the -algorithm at each step on www.simplexme.com (we used this to debug parts of our -code), but we found that SimplexMe only works for the ``simple'' case. It has -errors for large matrices and the infeasible and unbounded cases. - -\section{Planning} -Here is our \href{./DraftSpec_Annotated.pdf}{functional spec} and our -\href{./FinalSpec_Annotated.pdf}{technical spec} (the links open new PDFs). - -We planned to have a completed Simplex solver for simple cases, handling less -than or equal constraints in addition to a complete Matrix module. Here are our -initial functional spec and technical spec. We achieved far more than we -expected to achieve. - -Our matrix module is fairly robust and it certainly implements everything that -is necessary to program the simplex algorithm (which was our main goal for this -module). Yet, it goes further. We managed to implement determinant, trace, -row-reduction,inverse, LU decomposition, and more (all extra functions that we -classified under our ``cool features'' category. There were only minor points -that we did not achieve, such as the eigen function, which finds eigenvectors -and eigenvalues. - -As to the Elts module, we really went far and beyond. Because we coded -everything from an abstract point of view, once we had completed our project we -were able to go back and add a new Nums module which provided arbitrary -precision and size (ie, implementing arbitrary float arithmetic and bignum -arithmetic). Of course, to be honest, we really didn't actually code most of -this up. We used the built in Nums library. But as can be seen inside of -Elts.ml, our code was planned out so well that we can replace the foundation of -the mathematics in our program (\verb@ORDERED_AND_OPERATIONAL@) with only a few -changes in code. - -The Simplex module also went further than expected (though it certainly was more -difficult to implement than we expected). Not only did we implement the ability -to solve simple linear programs, but we can now solve with certainty any linear -program (dual and primal---dual requires \verb@initialize_simplex@ while primal -is the ``simple'' case), and provide the user with not only the optimal value, -but also one possible solution (which, in most cases is the only solution). - -Most of the milestones were achieved for all of our individual modules, with the -exception of our \verb@row_reduce@ (which was buggy and took longer than -expected), and of our simplex module (we went a day over what we had originally -planned). - -\section{Design and Implementation} - -Implementing simplex was difficult. A lot of difficulty stemmed from the -complexity and size of the algorithm itself, and we spent a lot of time trying -to understand each step. We never completely internalized every step of the -algorithm which made debugging difficult because it took a lot of effort to -recall the correct steps. We also had some initial terminology confused which -caused some of the code to be written incorrectly. Since each step ended up -being so complex, whenever we moved to the following steps of the algorithm we -always had trouble recalling what exactly we were doing and if someone else had -already done it. Furthermore, bugs were very difficult to find because the -system was so complex and there weren't many cases for which we could work -out the steps by hand. - -Finding SimplexMe was essential to finally getting all of the bugs out of our -system, but it also added many hours of unnecessary debugging. SimplexMe has -implemented simplex incorrectly and therefore the answer it gives for non-simple -systems can be, at times, incorrect. We verified this with Mathematica's -built in Linear Programming function. - -Testing the system also proved to be a problem in its own right. Even after we -got all of the parsing functions down and were certain that the data was being -loaded correctly, we were never really able to find a large set of simplex -example problems. In the end, we resorted to each team member individually -finding a problem and solution (with Mathematica, of course, since SimplexMe -proved unreliable), and then typing out the data into a text file so we could -test it. Once we had all the data, we decided to program an interface that would -allow for expanding the testability of the code. More information on this can be -found in the \verb@README@. - -To use the simplex algorithm, we first parse a standardized plain text file -containing the parameters of the objective and the constraints. The first line -of this file must either be \verb@max@ or \verb@min@, signalling a maximization -or minimization step. The next line denotes the objective function. The -coefficients of each of the variables in the objective are separated by columns -(decimals must be represented by fractions because we are using Ocaml's num -library). The last element in this is the ``constant'' (5 if the objective is -$2x + 3y + 5$). The next line must read \verb@subject to@ to help us with -parsing. The next lines are the constraints, where the coefficient of each -variable that appears in the objective must appear (in the same order as in the -objective), separated by commas. We also allow the user to enter if the -constraint is a $\geq, \leq$, or $=$ condition. Again, the last element is the -constant (5 in the case of $x + y \leq 5$). It is assumed that all the variables -in the objective function are $\geq 0$. - -Our program will parse this text file and generate a matrix representing the -objective and the constraints. Additionally, our program adds slack variables -(for example, (we change the inequalities in the constraints to equations---$3x -+ 2y \leq 5$ becomes $3x + 2y + s = 5$, where $s \geq 0$). This makes the -problem easier to solve). - -The simplex algorithm basically takes the objective function and attempts to -find a ``basic solution'' which will satisfy all the constraints. Even if this -solution is trivial (like setting all the variables in the objective to 0), we -take this solution and keep optimizing it by replacing a variable not in the -objective function with a variable that is. We check the objective to make sure -none of the coefficients are negative; if they are, we are finished and we have -the optimum value in the top right corner of the matrix (because we can set all -the variables in the objective to 0 and not violate the constraints). Our -``simple'' case is the case where we want to minimize the objective and all the -constraints are $\leq$ relations. However, we can change any problem to this -form. Maximization can be achieved by minimizing the negative of the objective. -Similarly, a $\geq$ constraint can be changed to a $\leq$ constraint by -multiplying each side of the inequality by $-1$. Finally, an $=$ constraint can -be changed into a $\geq$ and a $\leq$. - -As to the design and implementation of our actual code base, it definitely can -be improved. It is the best it can be given the fact that we hadn't already done -this same problem before, but as with anything, it can be better. For example, -our \verb@initialize_simplex@ function ended up being very complex and -complicated. We all understand the general idea of what it does, but not all of -us could implement it by ourselves from scratch. Basically, -\verb@initialize_simplex@ takes in a system and adds the slack variables to it. -During this step, if the algorithm determines a system (the objective -function and the constraints) to be unfeasible, it adds an extra variable and -attempts to minimize that (ie, replaces the objective function). If this new -system has a minimum of 0, then the original system is feasible and all we need -to do is substitute the old objective back in appropriately. Otherwise, the -original system is unfeasible and the program returns \verb@None@. In the first -case, an extra pivot operation might be necessary, but again, -\verb@initialize_simplex@ takes care of this. In the end, it calls -\verb@simple_solve@ (the function we originally wrote) on this now simple case -and a solution is produced. However, we did run our simplex algorithm against an -equation with 100 variables and 100 constraints (\verb@test22.txt@), and it -solved it in under 1 second. Even though the worst case running time of simplex -is not polynomial (it's exponential), on most real world cases simplex runs very -quickly, as evidenced by our test. - -Each member of the team contributed to the project. The work was split between -the matrix and elts modules, parsing the data, and the simplex algorithm. Andy -and Luis implemented most of the matrix module and parsing function, along with -the general framework of the program that allowed us to efficiently run and test -our code. Luis, Jason, and Ding worked on implementing and debugging the simplex -algorithm. Ding focused on comprehending every step of the algorithm as -presented in the Algorithms book and then presenting that information in a more -comprehensive manner to the group, while Jason focused on understanding the -process behind simplex and on coming up with example problems by hand. Jason and -Ding worked on testing simplex, and Luis and Andy provided much of the -functionality that allowed for easy testing. - -\section{Reflection} - -Our original planning went pretty well and we accomplished most of what we -expected to accomplish. We mostly hit both our milestones, though we missed our -second milestone by a day. We were definitely surprised by how difficult simplex -ended up being, and there were times when we thought we were essentially done -with the algorithm and it turned out that we were only halfway there. - -Given more time, we would like to polish up input and outputs of our algorithm, -including a more sophisticated parsing function and a user friendly output. -Within the matrix module itself, there are a few functions that we would have -liked to implement including finding the eigenvectors/eigenvalues and -calculating the norm of a matrix, among others. - -We made the right choice when we decided to drop some of our unnecessary matrix -module functions to spend more time on simplex, but we made an early mistake of -not bothering to understand simplex until the latter half of the project, and -rather than understanding the whole algorithm before we started coding, we chose -to only understand the steps we were directly implementing, leading to the -restructuring of major portions of our code towards the end. The next time we -tackle an algorithm-centric project like this, we should definitely plan out the -entire algorithm before we code (including all the possible edge cases that -might pop up). - -The most important thing we learned is how important good style, readability, -and careful planning is to programming and debugging. There were many times when -we were confused by what we wrote when we went back to debug our code, and since -no one was able to fully internalize every part of the algorithm, it was -important to abstract over helper functions and different parts of the code. -Having had someone who understood the algorithm entirely (or who had maybe -implemented something on this scale before), would have been extremely helpful -as it would have allowed us to better plan our time and resources as well as -abstract away many of the more detailed concepts. - -\section{Advice for Future Students} - -Start as early as possible and put some serious work in. Don't leave the video -until the last minute. This is a big project and showing it off is a very -rewarding part of it, especially when done well. Also try to understand the -general outline for every step before serious implementation, even if it means -delaying the implementation. Designing the framework is crucial and it makes the -whole process easier to get right the first time. Make sure to write test -functions as you write functions and to test each function right after you -finish writing it. This way, you will avoid having to hunt down a small bug in -one obscure part of your code. You'll know exactly when and where your code -went wrong. - -Additionally, make sure your group members don't play too much League of Legends -or Starcraft! - -But lastly---enjoy the project and be proud of what you've accomplished! - -\end{document} +\documentclass[letterpaper,12pt]{article} +% Change the header if you're going to change this. + +% Possible packages - uncomment to use +%\usepackage{amsmath} % Needed for math stuff +%\usepackage{lastpage} % Finds last +%\usepackage{amssymb} % Needed for some math symbols +%\usepackage{graphicx} % Needed for graphics +\usepackage[usenames,dvipsnames]{xcolor} % Needed for graphics and color +%\usepackage{setspace} % Needed to set single/double spacing +\usepackage{hyperref} + +%Suggested by TeXworks +\usepackage[utf8]{inputenc} % set input encoding (not needed with XeLaTeX) +\usepackage{verbatim} % adds environment for commenting out blocks of text & for better verbatim + +% Sets the page margins to 1 inch each side +\usepackage[margin=1in]{geometry} + +% Uncomment this section if you wish to have a header. +\usepackage{fancyhdr} +\pagestyle{fancy} +\renewcommand{\headrulewidth}{0.5pt} % customise the layout... +\lhead{CS51 Final Proect Report} \chead{} \rhead{May 5, 2013} +\lfoot{} \cfoot{\thepage} \rfoot{} + +\frenchspacing +\setlength\parindent{0pt} +\setlength\parskip{2ex} + +\newcommand{\annot}[1]{\textbf{\color{BrickRed} [#1]}} + +\newcommand{\vect}[1]{\mathbf{#1}} + +\begin{document} +\title{CS 51 Final Project Report} +\author{ +Luis Perez $|$ luisperez@college.harvard.edu \\ +Andy Shi $|$ andyshi@college.harvard.edu \\ +Zihao Wang $|$ zihaowang01@college.harvard.edu \\ +Ding Zhou $|$ dzhou@college.harvard.edu +} +\date{May 5, 2013} +\maketitle + +Link to demo video: \url{} + +\section{Overview} + +In the real world, people love to make money. Even in virtual worlds (bitcoin?), +people love to make money. Basically, many problems in business and even +everyday life (but especially business), boil down to maximizing profit or +minimizing costs. Of course, when you speak of minimization or maximization +there are always constraints. Resources are finite and time is precious. +Therefore, in an abstract sense, the problem boils down to minimizing or +maximizing some quantity under certain constraints. For example, you might want +to maximize the number of textbooks $x,y,z$ you buy under the constraints that you +can't buy more than five from one company, the total cost can't be more +than \$500, and the last two books you buy must cost more than \$75. This problem +is precisely the type of problem that our program can now solve. How did we do +this? + +To start off, we created a framework. We implemented a matrix module with +standard matrix operations such as multiplication, addition, map (where scaling +follows) as well as many extra features such as row reduction, determinant, +trace, and many more. The matrix module basically allows a user to do almost +anything imaginable with a matrix. In addition, we provided read and write +functionality for easy loading of data, and printing functionality for +straightforward viewing of the abstract matrix. However the intellectually +challenging part of our project is the implementation of the simplex algorithm. +We needed to create a matrix module simply because OCaml doesn't provide one +(also because we just wanted to implement a matrix module, of course). We use +many of the functions in the matrix module to create and manipulate simplex +systems within our simplex module. This allowed us to abstract away yet another +part of our project. + +The simplex algorithm is a type of linear programming that allows one to +optimize a linear equation given linear constraints (like the example problem at +the beginning). In mathematical terms, the simplex program operates on linear +programs in standard form. That is, it seeks to minimize (or maximize) +$C\vect{x}$ subject to $A\vect{x} = \vect{b} $ (Wikipedia). $\vect{x}$ and +$\vect{b}$ are vectors. $x$ contains the variables in the equation we want to +optimize (the objective), and $b$ contains the constants in the constraints. $C$ +and $A$ are matrices, and they represent the coefficients in the objective and +the constraints, respectively. Let's consider the inequalities: $30x + 15y + 18z +\geq 500, x \geq 5,y \geq 5, z \geq 5$, and $15y+18z \geq 75$ and the objective +functions $x +y + z$ where $x, y, z \geq 0$. As you will notice, this is +precisely the system of equations necessary to solve our book example above. The +algorithm converts the inequalities into equalities by adding positive slack +variables. Then it represents this set of equations in the matrix form where the +coefficients of the terms we maximize are on the first row and each column +represents a variable. This provides affordable storage of most of the +information necessary to solve the linear program. Furthermore, we then define +basic and non-basic variables (kept track of in a list), which together with the +aforementioned matrix, create a simplex system. + +The first function necessary for simplex is called \verb@simple_solve@. This +function, in abstract terms, finds the minimum value of a passed in system (as +described above). When we call the function \verb@simple_solve@, the program +finds the entering and leaving variables which are the variable that would +exchange position from basic and non-basic (slack to non-slack in the initial +case). Here, the algorithm is making the assumption that the system has a +feasible solution where all of our non-basic variables are zero, and therefore +exchanges the role of one of these variables with the corresponding basic +variable. In geometric terms, we are moving along one edge of the n-gon formed +by the constraints and attempting to minimize. This step is accomplished by the +helper function pivot. The function \verb@simple_solve@ continues to loop +(moving along edges and minimizing) until the first row---which represents the +term we want to optimize consists of only non-positive values or there is no +possible entering variable (ie, any change in the variables will only increase +our answer, not minimize it). In the first case, algebraically, it means that +the constraint can only increase if the variables change. Therefore, we have +found our minimum value which is the constant stored at the end of the first +row. This is the optimized value of the objective function. In the second case, +where no suitable entering variable can be found, it means that the system can +be decreased infinitely as pleased and therefore it is Unbounded. In the case +where a solution exists, one possible value for the original variables can also +be found by looking at their corresponding columns and checking to see if they +are basic. If there are an infinite number of these, our algorithm returns only +one. + +What if we had the constraint $x + y + z \geq -10$? As stated above, our +\verb@simple_solve@ makes the assumption that one of the feasible solutions is +the zero solution, but if we initially set $x, y$, and $z$ to 0 (our algorithm +generates a ``basic solution'' by doing this step) the constraint does not hold. +We therefore have a helper function called \verb@initialize_simplex@ that +converts everything into the ``simple'' form so we can just recursively call +\verb@simple_solve@. It does this by following a very complex algorithm better +detailed below. Our final simplex algorithm now accounts for all cases as long +as the equations are linear. It also handles not only less than or equal +constraints, but also greater than or equal to and equality constraints. +Finally, not only can it minimize, but it can also maximize. You can follow the +algorithm at each step on www.simplexme.com (we used this to debug parts of our +code), but we found that SimplexMe only works for the ``simple'' case. It has +errors for large matrices and the infeasible and unbounded cases. + +\section{Planning} +Here is our \href{./DraftSpec_Annotated.pdf}{functional spec} and our +\href{./FinalSpec_Annotated.pdf}{technical spec} (the links open new PDFs). + +We planned to have a completed Simplex solver for simple cases, handling less +than or equal constraints in addition to a complete Matrix module. Here are our +initial functional spec and technical spec. We achieved far more than we +expected to achieve. + +Our matrix module is fairly robust and it certainly implements everything that +is necessary to program the simplex algorithm (which was our main goal for this +module). Yet, it goes further. We managed to implement determinant, trace, +row-reduction,inverse, LU decomposition, and more (all extra functions that we +classified under our ``cool features'' category. There were only minor points +that we did not achieve, such as the eigen function, which finds eigenvectors +and eigenvalues. + +As to the Elts module, we really went far and beyond. Because we coded +everything from an abstract point of view, once we had completed our project we +were able to go back and add a new Nums module which provided arbitrary +precision and size (ie, implementing arbitrary float arithmetic and bignum +arithmetic). Of course, to be honest, we really didn't actually code most of +this up. We used the built in Nums library. But as can be seen inside of +Elts.ml, our code was planned out so well that we can replace the foundation of +the mathematics in our program (\verb@ORDERED_AND_OPERATIONAL@) with only a few +changes in code. + +The Simplex module also went further than expected (though it certainly was more +difficult to implement than we expected). Not only did we implement the ability +to solve simple linear programs, but we can now solve with certainty any linear +program (dual and primal---dual requires \verb@initialize_simplex@ while primal +is the ``simple'' case), and provide the user with not only the optimal value, +but also one possible solution (which, in most cases is the only solution). + +Most of the milestones were achieved for all of our individual modules, with the +exception of our \verb@row_reduce@ (which was buggy and took longer than +expected), and of our simplex module (we went a day over what we had originally +planned). + +\section{Design and Implementation} + +Implementing simplex was difficult. A lot of difficulty stemmed from the +complexity and size of the algorithm itself, and we spent a lot of time trying +to understand each step. We never completely internalized every step of the +algorithm which made debugging difficult because it took a lot of effort to +recall the correct steps. We also had some initial terminology confused which +caused some of the code to be written incorrectly. Since each step ended up +being so complex, whenever we moved to the following steps of the algorithm we +always had trouble recalling what exactly we were doing and if someone else had +already done it. Furthermore, bugs were very difficult to find because the +system was so complex and there weren't many cases for which we could work +out the steps by hand. + +Finding SimplexMe was essential to finally getting all of the bugs out of our +system, but it also added many hours of unnecessary debugging. SimplexMe has +implemented simplex incorrectly and therefore the answer it gives for non-simple +systems can be, at times, incorrect. We verified this with Mathematica's +built in Linear Programming function. + +Testing the system also proved to be a problem in its own right. Even after we +got all of the parsing functions down and were certain that the data was being +loaded correctly, we were never really able to find a large set of simplex +example problems. In the end, we resorted to each team member individually +finding a problem and solution (with Mathematica, of course, since SimplexMe +proved unreliable), and then typing out the data into a text file so we could +test it. Once we had all the data, we decided to program an interface that would +allow for expanding the testability of the code. More information on this can be +found in the \verb@README@. + +To use the simplex algorithm, we first parse a standardized plain text file +containing the parameters of the objective and the constraints. The first line +of this file must either be \verb@max@ or \verb@min@, signalling a maximization +or minimization step. The next line denotes the objective function. The +coefficients of each of the variables in the objective are separated by columns +(decimals must be represented by fractions because we are using Ocaml's num +library). The last element in this is the ``constant'' (5 if the objective is +$2x + 3y + 5$). The next line must read \verb@subject to@ to help us with +parsing. The next lines are the constraints, where the coefficient of each +variable that appears in the objective must appear (in the same order as in the +objective), separated by commas. We also allow the user to enter if the +constraint is a $\geq, \leq$, or $=$ condition. Again, the last element is the +constant (5 in the case of $x + y \leq 5$). It is assumed that all the variables +in the objective function are $\geq 0$. + +Our program will parse this text file and generate a matrix representing the +objective and the constraints. Additionally, our program adds slack variables +(for example, (we change the inequalities in the constraints to equations---$3x ++ 2y \leq 5$ becomes $3x + 2y + s = 5$, where $s \geq 0$). This makes the +problem easier to solve). + +The simplex algorithm basically takes the objective function and attempts to +find a ``basic solution'' which will satisfy all the constraints. Even if this +solution is trivial (like setting all the variables in the objective to 0), we +take this solution and keep optimizing it by replacing a variable not in the +objective function with a variable that is. We check the objective to make sure +none of the coefficients are negative; if they are, we are finished and we have +the optimum value in the top right corner of the matrix (because we can set all +the variables in the objective to 0 and not violate the constraints). Our +``simple'' case is the case where we want to minimize the objective and all the +constraints are $\leq$ relations. However, we can change any problem to this +form. Maximization can be achieved by minimizing the negative of the objective. +Similarly, a $\geq$ constraint can be changed to a $\leq$ constraint by +multiplying each side of the inequality by $-1$. Finally, an $=$ constraint can +be changed into a $\geq$ and a $\leq$. + +As to the design and implementation of our actual code base, it definitely can +be improved. It is the best it can be given the fact that we hadn't already done +this same problem before, but as with anything, it can be better. For example, +our \verb@initialize_simplex@ function ended up being very complex and +complicated. We all understand the general idea of what it does, but not all of +us could implement it by ourselves from scratch. Basically, +\verb@initialize_simplex@ takes in a system and adds the slack variables to it. +During this step, if the algorithm determines a system (the objective +function and the constraints) to be unfeasible, it adds an extra variable and +attempts to minimize that (ie, replaces the objective function). If this new +system has a minimum of 0, then the original system is feasible and all we need +to do is substitute the old objective back in appropriately. Otherwise, the +original system is unfeasible and the program returns \verb@None@. In the first +case, an extra pivot operation might be necessary, but again, +\verb@initialize_simplex@ takes care of this. In the end, it calls +\verb@simple_solve@ (the function we originally wrote) on this now simple case +and a solution is produced. However, we did run our simplex algorithm against an +equation with 100 variables and 100 constraints (\verb@test22.txt@), and it +solved it in under 1 second. Even though the worst case running time of simplex +is not polynomial (it's exponential), on most real world cases simplex runs very +quickly, as evidenced by our test. + +Each member of the team contributed to the project. The work was split between +the matrix and elts modules, parsing the data, and the simplex algorithm. Andy +and Luis implemented most of the matrix module and parsing function, along with +the general framework of the program that allowed us to efficiently run and test +our code. Luis, Jason, and Ding worked on implementing and debugging the simplex +algorithm. Ding focused on comprehending every step of the algorithm as +presented in the Algorithms book and then presenting that information in a more +comprehensive manner to the group, while Jason focused on understanding the +process behind simplex and on coming up with example problems by hand. Jason and +Ding worked on testing simplex, and Luis and Andy provided much of the +functionality that allowed for easy testing. + +\section{Reflection} + +Our original planning went pretty well and we accomplished most of what we +expected to accomplish. We mostly hit both our milestones, though we missed our +second milestone by a day. We were definitely surprised by how difficult simplex +ended up being, and there were times when we thought we were essentially done +with the algorithm and it turned out that we were only halfway there. + +Given more time, we would like to polish up input and outputs of our algorithm, +including a more sophisticated parsing function and a user friendly output. +Within the matrix module itself, there are a few functions that we would have +liked to implement including finding the eigenvectors/eigenvalues and +calculating the norm of a matrix, among others. + +We made the right choice when we decided to drop some of our unnecessary matrix +module functions to spend more time on simplex, but we made an early mistake of +not bothering to understand simplex until the latter half of the project, and +rather than understanding the whole algorithm before we started coding, we chose +to only understand the steps we were directly implementing, leading to the +restructuring of major portions of our code towards the end. The next time we +tackle an algorithm-centric project like this, we should definitely plan out the +entire algorithm before we code (including all the possible edge cases that +might pop up). + +The most important thing we learned is how important good style, readability, +and careful planning is to programming and debugging. There were many times when +we were confused by what we wrote when we went back to debug our code, and since +no one was able to fully internalize every part of the algorithm, it was +important to abstract over helper functions and different parts of the code. +Having had someone who understood the algorithm entirely (or who had maybe +implemented something on this scale before), would have been extremely helpful +as it would have allowed us to better plan our time and resources as well as +abstract away many of the more detailed concepts. + +\section{Advice for Future Students} + +Start as early as possible and put some serious work in. Don't leave the video +until the last minute. This is a big project and showing it off is a very +rewarding part of it, especially when done well. Also try to understand the +general outline for every step before serious implementation, even if it means +delaying the implementation. Designing the framework is crucial and it makes the +whole process easier to get right the first time. Make sure to write test +functions as you write functions and to test each function right after you +finish writing it. This way, you will avoid having to hunt down a small bug in +one obscure part of your code. You'll know exactly when and where your code +went wrong. + +Additionally, make sure your group members don't play too much League of Legends +or Starcraft! + +But lastly---enjoy the project and be proud of what you've accomplished! + +\end{document} diff --git a/tests/simplex/results.txt b/tests/simplex/results.txt deleted file mode 100644 index f0855f0..0000000 --- a/tests/simplex/results.txt +++ /dev/null @@ -1,22 +0,0 @@ -None -Unbounded -28,(0, 1/6, 2/3) -22,(3, 0, 0) -0,(0, 0, 0) -Unbounded -3,(0, 3/2) -2,(14/9, 10/9) -5,(0, 0, 0) -65,(0, 69/7, 2/7, 0) --130/7,(15/7, 0, 25/7) -33,(18, 5) --39901/3856,(0, 0, 0, 4137/1928, 197/964, 0, 4623/3856, 103/3856) -Unbounded -1000000,(20, 20, 0) -960000,(0, 40, 0) -2500,(160/11, 210/11, 0) -2650,(0, 0, 24) -20100,(140, 130) -1080,(72, 0, 0) -33,(3,12) -1,(1, 0) \ No newline at end of file diff --git a/tests/simplex/test0.txt b/tests/simplex/test0.txt deleted file mode 100644 index 1ccbe3c..0000000 --- a/tests/simplex/test0.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -1,-2, 0 -subject to -1,2,<=,4 --2,-6,<=,-12 -0,1,<=,1 diff --git a/tests/simplex/test1.txt b/tests/simplex/test1.txt deleted file mode 100644 index 6a6bb4f..0000000 --- a/tests/simplex/test1.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -1,3,0 -subject to --1,1,<=,-1 --1,-1,<=,-3 --1,4,<=,2 \ No newline at end of file diff --git a/tests/simplex/test10.txt b/tests/simplex/test10.txt deleted file mode 100644 index 34a4562..0000000 --- a/tests/simplex/test10.txt +++ /dev/null @@ -1,5 +0,0 @@ -min --2,-3,-4,0 -subject to -3,2,1,=,10 -2,5,3,=,15 \ No newline at end of file diff --git a/tests/simplex/test11.txt b/tests/simplex/test11.txt deleted file mode 100644 index ebd8c3f..0000000 --- a/tests/simplex/test11.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -1,3,0 -subject to -1,-2,<=,8 --1,-1,<=,-3 --1,4,<=,2 \ No newline at end of file diff --git a/tests/simplex/test12.txt b/tests/simplex/test12.txt deleted file mode 100644 index 41a72fc..0000000 --- a/tests/simplex/test12.txt +++ /dev/null @@ -1,11 +0,0 @@ -min -5,7,23,-5,-10,7,2,1,0 -subject to -0,0,0,2,-5,-1,-6,7,<=,1 -21,5,5,0,15,1,-1,5,<=,2 -2,3,5,-5,-6,-1,-1,-5,<=,3 --4,-5,-3,-6,7,-4,1,3,<=,4 -0,0,2,5,3,-1,-5,-13,<=,5 -0,1,2,4,3,1,-3,15,<=,6 -5,1,3,2,-5,1,3,5,<=,7 -14,13,52,13,-41,12,14,-31,<=,50 \ No newline at end of file diff --git a/tests/simplex/test13.txt b/tests/simplex/test13.txt deleted file mode 100644 index c6988b7..0000000 --- a/tests/simplex/test13.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -6,8,0 -subject to --1,2,<=,6 --1,1,<=,7 -0,1,<=,8 \ No newline at end of file diff --git a/tests/simplex/test14.txt b/tests/simplex/test14.txt deleted file mode 100644 index 8790df8..0000000 --- a/tests/simplex/test14.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -8000,6000,-24000,720000 -subject to -1/2,1/2,1,<=,30 -2,2,-8,<=,80 -2,1,-4,<=,60 \ No newline at end of file diff --git a/tests/simplex/test15.txt b/tests/simplex/test15.txt deleted file mode 100644 index b98b800..0000000 --- a/tests/simplex/test15.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -5000,6000,-24000,720000 -subject to -1/2,1/2,1,<=,30 -2,2,-8,<=,80 -2,1,-4,<=,60 \ No newline at end of file diff --git a/tests/simplex/test16.txt b/tests/simplex/test16.txt deleted file mode 100644 index 2c3e0f7..0000000 --- a/tests/simplex/test16.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -80,70,100,0 -subject to -3,4,5,<=,120 -5,3,5,<=,130 -4,3,5,<=,120 \ No newline at end of file diff --git a/tests/simplex/test17.txt b/tests/simplex/test17.txt deleted file mode 100644 index 56e47d3..0000000 --- a/tests/simplex/test17.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -80,70,110,0 -subject to -3,4,5,<=,120 -5,3,5,<=,130 -4,3,5,<=,120 \ No newline at end of file diff --git a/tests/simplex/test18.txt b/tests/simplex/test18.txt deleted file mode 100644 index d574658..0000000 --- a/tests/simplex/test18.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -60,90,0 -subject to -1,1,<=,320 -50,100,<=,20000 -100,40,<=,19200 \ No newline at end of file diff --git a/tests/simplex/test19.txt b/tests/simplex/test19.txt deleted file mode 100644 index b0b60ba..0000000 --- a/tests/simplex/test19.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -15,18,20,0 -subject to -1/6,1/4,1/2,<=,60 -40,50,60,<=,2880 -25,30,40,<=,2400 \ No newline at end of file diff --git a/tests/simplex/test2.txt b/tests/simplex/test2.txt deleted file mode 100644 index c49edb0..0000000 --- a/tests/simplex/test2.txt +++ /dev/null @@ -1,6 +0,0 @@ -min -30,24,36,0 -subject to -1,2,4,>=,3 -1,2,1,>=,1 -3,5,2,>=,2 \ No newline at end of file diff --git a/tests/simplex/test20.txt b/tests/simplex/test20.txt deleted file mode 100644 index fce80e8..0000000 --- a/tests/simplex/test20.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -3,2,0 -subject to -2,1,<=,18 -2,3,<=,42 -3,1,<=,24 \ No newline at end of file diff --git a/tests/simplex/test21.txt b/tests/simplex/test21.txt deleted file mode 100644 index 42d4372..0000000 --- a/tests/simplex/test21.txt +++ /dev/null @@ -1,4 +0,0 @@ -max -1,1,0 -subject to -1,1,<=,1 \ No newline at end of file diff --git a/tests/simplex/test3.txt b/tests/simplex/test3.txt deleted file mode 100644 index c4b62c4..0000000 --- a/tests/simplex/test3.txt +++ /dev/null @@ -1,5 +0,0 @@ -max -4,3,2,10 -subject to -1,2,3,<=,3 -4,5,6,<=,14 \ No newline at end of file diff --git a/tests/simplex/test4.txt b/tests/simplex/test4.txt deleted file mode 100644 index f30621c..0000000 --- a/tests/simplex/test4.txt +++ /dev/null @@ -1,5 +0,0 @@ -min -1,1,1,0 -subject to -1,1,0,<=,8 -0,-1,1,<=,0 \ No newline at end of file diff --git a/tests/simplex/test5.txt b/tests/simplex/test5.txt deleted file mode 100644 index e4e67f4..0000000 --- a/tests/simplex/test5.txt +++ /dev/null @@ -1,4 +0,0 @@ -max -1,1,0 -subject to -1,-2,<=,-1 \ No newline at end of file diff --git a/tests/simplex/test6.txt b/tests/simplex/test6.txt deleted file mode 100644 index 784a5ed..0000000 --- a/tests/simplex/test6.txt +++ /dev/null @@ -1,5 +0,0 @@ -min -1,2,0 -subject to -1,-2,<=,-3 -4,1,<=,2 \ No newline at end of file diff --git a/tests/simplex/test7.txt b/tests/simplex/test7.txt deleted file mode 100644 index c300d79..0000000 --- a/tests/simplex/test7.txt +++ /dev/null @@ -1,5 +0,0 @@ -max -2,-1,0 -subject to -2,-1,<=,2 -1,-5,<=,-4 \ No newline at end of file diff --git a/tests/simplex/test8.txt b/tests/simplex/test8.txt deleted file mode 100644 index 2ba0eb9..0000000 --- a/tests/simplex/test8.txt +++ /dev/null @@ -1,5 +0,0 @@ -max --2,-3,-4,5 -subject to -3,2,1,<=,10 -2,5,3,<=,15 \ No newline at end of file diff --git a/tests/simplex/test9.txt b/tests/simplex/test9.txt deleted file mode 100644 index a52c30f..0000000 --- a/tests/simplex/test9.txt +++ /dev/null @@ -1,6 +0,0 @@ -max -5,6,-4,1,7 -subject to -1,1,0,0,>=,-6 -0,0,7,-1,=,2 -2,2,1,1,<=,20 \ No newline at end of file