1:   2:   3:   4:   5:   6:   7:   8:   9:  10:  11:  12:  13:  14:  15:  16:  17:  18:  19:  20:  21:  22:  23:  24:  25:  26:  27:  28:  29:  30:  31:  32:  33:  34:  35:  36:  37:  38:  39:  40:  41:  42:  43:  44:  45:  46:  47:  48:  49:  50:  51:  52:  53:  54:  55:  56:  57:  58:  59:  60:  61:  62:  63:  64:  65:  66:  67:  68:  69:  70:  71:  72:  73:  74:  75:  76:  77:  78:  79:  80:  81:  82:  83:  84:  85:  86:  87:  88:  89:  90:  91:  92:  93:  94:  95:  96:  97:  98:  99: 100: 101: 102: 103: 104: 105: 106: 107: 108: 109: 110: 111: 112: 113: 114: 115: 116: 117: 118: 119: 120: 121: 122: 123: 124: 125: 126: 127: 128: 129: 130: 131: 132: 133: 134: 135: 136: 137: 138: 139: 140: 141: 142: 143: 144: 145: 146: 147: 148: 149: 150: 151: 152: 153: 154: 155: 156: 157: 158: 159: 160: 161: 162: 163: 164: 165: 166: 167: 168: 169: 170: 171: 172: 173: 174: 175: 176: 177: 178: 179: 180: 181: 182: 183: 184: 185: 186: 187: 188: 189: 190: 191: 192: 193: 194: 195: 196: 197: 198: 199: 200: 201: 202: 203: 204: 205: 206: 207: 208: 209: 210: 211: 212: 213: 214: 215: 216: 217: 218: 219: 220: 221: 222: 223: 224: 

<?php
/*
 * Math-Matrix library
 *
 * @author Ashley Kitson <akitson@zf4.biz>
 * @copyright Ashley Kitson, UK, 2014
 * @licence GPL V3 or later : http://www.gnu.org/licenses/gpl.html
 * @link http://en.wikipedia.org/wiki/Matrix_(mathematics)#Determinant
 * @link http://en.wikipedia.org/wiki/Laplace_expansion
 */

namespace Chippyash\Math\Matrix\Derivative\Strategy\Determinant;

use Chippyash\Math\Matrix\NumericMatrix;
use Chippyash\Math\Matrix\Interfaces\DeterminantStrategyInterface;
use Chippyash\Matrix\Transformation\Cofactor;
use Chippyash\Math\Type\Calculator;
use Chippyash\Math\Matrix\Traits\CreateCorrectScalarType;
use Chippyash\Math\Matrix\Interfaces\TuningInterface;

/**
 * Laplace expansion strategy for matrix determinant
 * Computing up to about M[9] is about as much as you can
 * do with this method of computing determinants in a reasonable time.
 *
 * @link examples/example-laplace-bounds.php
 *
 */
class Laplace implements DeterminantStrategyInterface, TuningInterface
{
    use CreateCorrectScalarType;

    /**
     * Cofactor function
     * @var Chippyash\Matrix\Transformation\Cofactor
     */
    protected $fCof;

    /**
     * Calculator
     * @var Chippyash\Math\Type\Calculator
     */
    protected $calc;

    /**
     * Cache for part processed determinants
     *
     * @var array
     */
    protected static $cache = [];

    /**
     * Compute determinant using laplace expansion
     * $mA must be
     * - square
     * This is not checked here - that is done in the determinant derivative class
     *
     * @param \Chippyash\Math\Matrix\NumericMatrix $mA
     * @return \Chippyash\Type\Interfaces\NumericTypeInterface
     */
    public function determinant(NumericMatrix $mA)
    {
        return $this->det($mA);
    }

    /**
     * Tune an item on a class. Available items:
     * - clearCache : boolean - if $value == true, will clear determinant cache
     *                Always returns number of items currently in cache
     *
     * @param \Chippyash\Type\String\StringType $name Item to tune
     * @param mixed $value Value to set
     *
     * @return mixed - previous value of item
     *
     * @throws \InvalidArgumentException if name does not exist
     */
    public function tune(\Chippyash\Type\String\StringType $name, $value)
    {
        if ($name() !== 'clearCache') {
            throw new \InvalidArgumentException("{$name} is unknown for tuning");
        }

        $ret = count(self::$cache);
        if ($value) {
            self::$cache = [];
        }
        return $ret;
    }

    /**
     * Recursive determinant function
     *
     * @param \Chippyash\Math\Matrix\NumericMatrix $mA
     * @return \Chippyash\Type\Interfaces\NumericTypeInterface
     */
    protected function det(NumericMatrix $mA)
    {
        $rowCount = $mA->rows();

        if ($rowCount == 0) {
            return $this->createCorrectScalarType($mA, 1);
        }

        if ($rowCount == 1) {
            return $mA->get(1,1);
        }

        $possAnswer = $this->checkInCache($mA);
        if ($possAnswer !== false) {
            return $possAnswer;
        }

        if ($rowCount == 2) {
            //2X2 matrix
            $det = $this->det2($mA);
        } else {
            //nXn matrix
            $det = $this->detN($mA);
        }
        $this->storeInCache($mA, $det);

        return $det;
    }

    /**
     * Return determinant of a 2X2 matrix
     * @link http://en.wikipedia.org/wiki/Matrix_determinant#2.C2.A0.C3.97.C2.A02_matrices
     * [[a, b]
     *  [c, d]]
     * ad - bc
     *
     * @param \Chippyash\Math\Matrix\NumericMatrix $mA
     */
    protected function det2(NumericMatrix $mA)
    {
        $c = $this->calc();
        return $c->sub(
                    $c->mul($mA->get(1,1), $mA->get(2,2)),
                    $c->mul($mA->get(1,2), $mA->get(2,1))
                );
    }

    /**
     * Get determinant for arbitrarily large matrices
     *
     * @link http://www.intmath.com/matrices-determinants/2-large-determinants.php
     * @param \Chippyash\Math\Matrix\NumericMatrix $mA
     * @return Chippyash\Type\Interfaces\NumericTypeInterface
     */
    protected function detN(NumericMatrix $mA)
    {

        $det = $this->createCorrectScalarType($mA, 0);
        $positive = true;
        $rowLimit = $mA->rows() + 1;
        $c = $this->calc();
        for ($r = 1; $r < $rowLimit; $r ++) {
            $t = $c->mul(
                    $mA->get($r, 1),
                    $this->det($this->cof()->transform($mA,array($r, 1)))
                    );
            $det = $c->add(
                    $det, ($positive ? $t : $t->negate())
                    );
            $positive = !$positive;
        }

        return $det;
    }

    /**
     * @return Chippyash\Matrix\Transformation\Cofactor
     */
    protected function cof()
    {
        if (empty($this->fCof)) {
            $this->fCof = new Cofactor();
        }

        return $this->fCof;
    }

    /**
     * @return Chippyash\Math\Type\Calculator
     */
    protected function calc()
    {
        if (empty($this->calc)) {
            $this->calc = new Calculator();
        }

        return $this->calc;
    }

    /**
     * Check for matrix determinant answer in cache
     *
     * @param NumericMatrix $mA
     * @return boolean|numeric
     */
    protected function checkInCache($mA)
    {
        $key = md5(serialize($mA));
        if (isset(self::$cache[$key])) {
            return self::$cache[$key];
        }

        return false;
    }

    /**
     * Store matrix determinant answer in cache
     *
     * @param NumericMatrix $mA
     * @param numeric $answer
     */
    protected function storeInCache($mA, $answer)
    {
        $key = md5(serialize($mA));
        self::$cache[$key] = $answer;
    }
}