MDL-58859 mlbackend_php: Added to core

Part of MDL-57791 epic.
This commit is contained in:
David Monllao
2017-07-24 07:51:54 +02:00
parent 229ae617c5
commit 9e84757bf0
129 changed files with 11485 additions and 0 deletions
@@ -0,0 +1,345 @@
<?php
declare(strict_types=1);
namespace Phpml\Association;
use Phpml\Helper\Predictable;
use Phpml\Helper\Trainable;
class Apriori implements Associator
{
use Trainable, Predictable;
const ARRAY_KEY_ANTECEDENT = 'antecedent';
const ARRAY_KEY_CONFIDENCE = 'confidence';
const ARRAY_KEY_CONSEQUENT = 'consequent';
const ARRAY_KEY_SUPPORT = 'support';
/**
* Minimum relative probability of frequent transactions.
*
* @var float
*/
private $confidence;
/**
* The large set contains frequent k-length item sets.
*
* @var mixed[][][]
*/
private $large;
/**
* Minimum relative frequency of transactions.
*
* @var float
*/
private $support;
/**
* The generated Apriori association rules.
*
* @var mixed[][]
*/
private $rules;
/**
* Apriori constructor.
*
* @param float $support
* @param float $confidence
*/
public function __construct(float $support = 0.0, float $confidence = 0.0)
{
$this->support = $support;
$this->confidence = $confidence;
}
/**
* Get all association rules which are generated for every k-length frequent item set.
*
* @return mixed[][]
*/
public function getRules() : array
{
if (!$this->large) {
$this->large = $this->apriori();
}
if ($this->rules) {
return $this->rules;
}
$this->rules = [];
$this->generateAllRules();
return $this->rules;
}
/**
* Generates frequent item sets.
*
* @return mixed[][][]
*/
public function apriori() : array
{
$L = [];
$L[1] = $this->items();
$L[1] = $this->frequent($L[1]);
for ($k = 2; !empty($L[$k - 1]); ++$k) {
$L[$k] = $this->candidates($L[$k - 1]);
$L[$k] = $this->frequent($L[$k]);
}
return $L;
}
/**
* @param mixed[] $sample
*
* @return mixed[][]
*/
protected function predictSample(array $sample) : array
{
$predicts = array_values(array_filter($this->getRules(), function ($rule) use ($sample) {
return $this->equals($rule[self::ARRAY_KEY_ANTECEDENT], $sample);
}));
return array_map(function ($rule) {
return $rule[self::ARRAY_KEY_CONSEQUENT];
}, $predicts);
}
/**
* Generate rules for each k-length frequent item set.
*/
private function generateAllRules()
{
for ($k = 2; !empty($this->large[$k]); ++$k) {
foreach ($this->large[$k] as $frequent) {
$this->generateRules($frequent);
}
}
}
/**
* Generate confident rules for frequent item set.
*
* @param mixed[] $frequent
*/
private function generateRules(array $frequent)
{
foreach ($this->antecedents($frequent) as $antecedent) {
if ($this->confidence <= ($confidence = $this->confidence($frequent, $antecedent))) {
$consequent = array_values(array_diff($frequent, $antecedent));
$this->rules[] = [
self::ARRAY_KEY_ANTECEDENT => $antecedent,
self::ARRAY_KEY_CONSEQUENT => $consequent,
self::ARRAY_KEY_SUPPORT => $this->support($consequent),
self::ARRAY_KEY_CONFIDENCE => $confidence,
];
}
}
}
/**
* Generates the power set for given item set $sample.
*
* @param mixed[] $sample
*
* @return mixed[][]
*/
private function powerSet(array $sample) : array
{
$results = [[]];
foreach ($sample as $item) {
foreach ($results as $combination) {
$results[] = array_merge([$item], $combination);
}
}
return $results;
}
/**
* Generates all proper subsets for given set $sample without the empty set.
*
* @param mixed[] $sample
*
* @return mixed[][]
*/
private function antecedents(array $sample) : array
{
$cardinality = count($sample);
$antecedents = $this->powerSet($sample);
return array_filter($antecedents, function ($antecedent) use ($cardinality) {
return (count($antecedent) != $cardinality) && ($antecedent != []);
});
}
/**
* Calculates frequent k = 1 item sets.
*
* @return mixed[][]
*/
private function items() : array
{
$items = [];
foreach ($this->samples as $sample) {
foreach ($sample as $item) {
if (!in_array($item, $items, true)) {
$items[] = $item;
}
}
}
return array_map(function ($entry) {
return [$entry];
}, $items);
}
/**
* Returns frequent item sets only.
*
* @param mixed[][] $samples
*
* @return mixed[][]
*/
private function frequent(array $samples) : array
{
return array_filter($samples, function ($entry) {
return $this->support($entry) >= $this->support;
});
}
/**
* Calculates frequent k item sets, where count($samples) == $k - 1.
*
* @param mixed[][] $samples
*
* @return mixed[][]
*/
private function candidates(array $samples) : array
{
$candidates = [];
foreach ($samples as $p) {
foreach ($samples as $q) {
if (count(array_merge(array_diff($p, $q), array_diff($q, $p))) != 2) {
continue;
}
$candidate = array_unique(array_merge($p, $q));
if ($this->contains($candidates, $candidate)) {
continue;
}
foreach ((array) $this->samples as $sample) {
if ($this->subset($sample, $candidate)) {
$candidates[] = $candidate;
continue 2;
}
}
}
}
return $candidates;
}
/**
* Calculates confidence for $set. Confidence is the relative amount of sets containing $subset which also contain
* $set.
*
* @param mixed[] $set
* @param mixed[] $subset
*
* @return float
*/
private function confidence(array $set, array $subset) : float
{
return $this->support($set) / $this->support($subset);
}
/**
* Calculates support for item set $sample. Support is the relative amount of sets containing $sample in the data
* pool.
*
* @see \Phpml\Association\Apriori::samples
*
* @param mixed[] $sample
*
* @return float
*/
private function support(array $sample) : float
{
return $this->frequency($sample) / count($this->samples);
}
/**
* Counts occurrences of $sample as subset in data pool.
*
* @see \Phpml\Association\Apriori::samples
*
* @param mixed[] $sample
*
* @return int
*/
private function frequency(array $sample) : int
{
return count(array_filter($this->samples, function ($entry) use ($sample) {
return $this->subset($entry, $sample);
}));
}
/**
* Returns true if set is an element of system.
*
* @see \Phpml\Association\Apriori::equals()
*
* @param mixed[][] $system
* @param mixed[] $set
*
* @return bool
*/
private function contains(array $system, array $set) : bool
{
return (bool) array_filter($system, function ($entry) use ($set) {
return $this->equals($entry, $set);
});
}
/**
* Returns true if subset is a (proper) subset of set by its items string representation.
*
* @param mixed[] $set
* @param mixed[] $subset
*
* @return bool
*/
private function subset(array $set, array $subset) : bool
{
return !array_diff($subset, array_intersect($subset, $set));
}
/**
* Returns true if string representation of items does not differ.
*
* @param mixed[] $set1
* @param mixed[] $set2
*
* @return bool
*/
private function equals(array $set1, array $set2) : bool
{
return array_diff($set1, $set2) == array_diff($set2, $set1);
}
}
@@ -0,0 +1,11 @@
<?php
declare(strict_types=1);
namespace Phpml\Association;
use Phpml\Estimator;
interface Associator extends Estimator
{
}
@@ -0,0 +1,11 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification;
use Phpml\Estimator;
interface Classifier extends Estimator
{
}
@@ -0,0 +1,509 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification;
use Phpml\Exception\InvalidArgumentException;
use Phpml\Helper\Predictable;
use Phpml\Helper\Trainable;
use Phpml\Math\Statistic\Mean;
use Phpml\Classification\DecisionTree\DecisionTreeLeaf;
class DecisionTree implements Classifier
{
use Trainable, Predictable;
const CONTINUOUS = 1;
const NOMINAL = 2;
/**
* @var array
*/
protected $columnTypes;
/**
* @var array
*/
private $labels = [];
/**
* @var int
*/
private $featureCount = 0;
/**
* @var DecisionTreeLeaf
*/
protected $tree = null;
/**
* @var int
*/
protected $maxDepth;
/**
* @var int
*/
public $actualDepth = 0;
/**
* @var int
*/
private $numUsableFeatures = 0;
/**
* @var array
*/
private $selectedFeatures;
/**
* @var array
*/
private $featureImportances = null;
/**
*
* @var array
*/
private $columnNames = null;
/**
* @param int $maxDepth
*/
public function __construct(int $maxDepth = 10)
{
$this->maxDepth = $maxDepth;
}
/**
* @param array $samples
* @param array $targets
*/
public function train(array $samples, array $targets)
{
$this->samples = array_merge($this->samples, $samples);
$this->targets = array_merge($this->targets, $targets);
$this->featureCount = count($this->samples[0]);
$this->columnTypes = self::getColumnTypes($this->samples);
$this->labels = array_keys(array_count_values($this->targets));
$this->tree = $this->getSplitLeaf(range(0, count($this->samples) - 1));
// Each time the tree is trained, feature importances are reset so that
// we will have to compute it again depending on the new data
$this->featureImportances = null;
// If column names are given or computed before, then there is no
// need to init it and accidentally remove the previous given names
if ($this->columnNames === null) {
$this->columnNames = range(0, $this->featureCount - 1);
} elseif (count($this->columnNames) > $this->featureCount) {
$this->columnNames = array_slice($this->columnNames, 0, $this->featureCount);
} elseif (count($this->columnNames) < $this->featureCount) {
$this->columnNames = array_merge($this->columnNames,
range(count($this->columnNames), $this->featureCount - 1));
}
}
/**
* @param array $samples
* @return array
*/
public static function getColumnTypes(array $samples) : array
{
$types = [];
$featureCount = count($samples[0]);
for ($i=0; $i < $featureCount; $i++) {
$values = array_column($samples, $i);
$isCategorical = self::isCategoricalColumn($values);
$types[] = $isCategorical ? self::NOMINAL : self::CONTINUOUS;
}
return $types;
}
/**
* @param array $records
* @param int $depth
* @return DecisionTreeLeaf
*/
protected function getSplitLeaf(array $records, int $depth = 0) : DecisionTreeLeaf
{
$split = $this->getBestSplit($records);
$split->level = $depth;
if ($this->actualDepth < $depth) {
$this->actualDepth = $depth;
}
// Traverse all records to see if all records belong to the same class,
// otherwise group the records so that we can classify the leaf
// in case maximum depth is reached
$leftRecords = [];
$rightRecords= [];
$remainingTargets = [];
$prevRecord = null;
$allSame = true;
foreach ($records as $recordNo) {
// Check if the previous record is the same with the current one
$record = $this->samples[$recordNo];
if ($prevRecord && $prevRecord != $record) {
$allSame = false;
}
$prevRecord = $record;
// According to the split criteron, this record will
// belong to either left or the right side in the next split
if ($split->evaluate($record)) {
$leftRecords[] = $recordNo;
} else {
$rightRecords[]= $recordNo;
}
// Group remaining targets
$target = $this->targets[$recordNo];
if (! array_key_exists($target, $remainingTargets)) {
$remainingTargets[$target] = 1;
} else {
$remainingTargets[$target]++;
}
}
if ($allSame || $depth >= $this->maxDepth || count($remainingTargets) === 1) {
$split->isTerminal = 1;
arsort($remainingTargets);
$split->classValue = key($remainingTargets);
} else {
if ($leftRecords) {
$split->leftLeaf = $this->getSplitLeaf($leftRecords, $depth + 1);
}
if ($rightRecords) {
$split->rightLeaf= $this->getSplitLeaf($rightRecords, $depth + 1);
}
}
return $split;
}
/**
* @param array $records
* @return DecisionTreeLeaf
*/
protected function getBestSplit(array $records) : DecisionTreeLeaf
{
$targets = array_intersect_key($this->targets, array_flip($records));
$samples = array_intersect_key($this->samples, array_flip($records));
$samples = array_combine($records, $this->preprocess($samples));
$bestGiniVal = 1;
$bestSplit = null;
$features = $this->getSelectedFeatures();
foreach ($features as $i) {
$colValues = [];
foreach ($samples as $index => $row) {
$colValues[$index] = $row[$i];
}
$counts = array_count_values($colValues);
arsort($counts);
$baseValue = key($counts);
$gini = $this->getGiniIndex($baseValue, $colValues, $targets);
if ($bestSplit === null || $bestGiniVal > $gini) {
$split = new DecisionTreeLeaf();
$split->value = $baseValue;
$split->giniIndex = $gini;
$split->columnIndex = $i;
$split->isContinuous = $this->columnTypes[$i] == self::CONTINUOUS;
$split->records = $records;
// If a numeric column is to be selected, then
// the original numeric value and the selected operator
// will also be saved into the leaf for future access
if ($this->columnTypes[$i] == self::CONTINUOUS) {
$matches = [];
preg_match("/^([<>=]{1,2})\s*(.*)/", strval($split->value), $matches);
$split->operator = $matches[1];
$split->numericValue = floatval($matches[2]);
}
$bestSplit = $split;
$bestGiniVal = $gini;
}
}
return $bestSplit;
}
/**
* Returns available features/columns to the tree for the decision making
* process. <br>
*
* If a number is given with setNumFeatures() method, then a random selection
* of features up to this number is returned. <br>
*
* If some features are manually selected by use of setSelectedFeatures(),
* then only these features are returned <br>
*
* If any of above methods were not called beforehand, then all features
* are returned by default.
*
* @return array
*/
protected function getSelectedFeatures() : array
{
$allFeatures = range(0, $this->featureCount - 1);
if ($this->numUsableFeatures === 0 && ! $this->selectedFeatures) {
return $allFeatures;
}
if ($this->selectedFeatures) {
return $this->selectedFeatures;
}
$numFeatures = $this->numUsableFeatures;
if ($numFeatures > $this->featureCount) {
$numFeatures = $this->featureCount;
}
shuffle($allFeatures);
$selectedFeatures = array_slice($allFeatures, 0, $numFeatures, false);
sort($selectedFeatures);
return $selectedFeatures;
}
/**
* @param $baseValue
* @param array $colValues
* @param array $targets
* @return float
*/
public function getGiniIndex($baseValue, array $colValues, array $targets) : float
{
$countMatrix = [];
foreach ($this->labels as $label) {
$countMatrix[$label] = [0, 0];
}
foreach ($colValues as $index => $value) {
$label = $targets[$index];
$rowIndex = $value === $baseValue ? 0 : 1;
$countMatrix[$label][$rowIndex]++;
}
$giniParts = [0, 0];
for ($i=0; $i<=1; $i++) {
$part = 0;
$sum = array_sum(array_column($countMatrix, $i));
if ($sum > 0) {
foreach ($this->labels as $label) {
$part += pow($countMatrix[$label][$i] / floatval($sum), 2);
}
}
$giniParts[$i] = (1 - $part) * $sum;
}
return array_sum($giniParts) / count($colValues);
}
/**
* @param array $samples
* @return array
*/
protected function preprocess(array $samples) : array
{
// Detect and convert continuous data column values into
// discrete values by using the median as a threshold value
$columns = [];
for ($i=0; $i<$this->featureCount; $i++) {
$values = array_column($samples, $i);
if ($this->columnTypes[$i] == self::CONTINUOUS) {
$median = Mean::median($values);
foreach ($values as &$value) {
if ($value <= $median) {
$value = "<= $median";
} else {
$value = "> $median";
}
}
}
$columns[] = $values;
}
// Below method is a strange yet very simple & efficient method
// to get the transpose of a 2D array
return array_map(null, ...$columns);
}
/**
* @param array $columnValues
* @return bool
*/
protected static function isCategoricalColumn(array $columnValues) : bool
{
$count = count($columnValues);
// There are two main indicators that *may* show whether a
// column is composed of discrete set of values:
// 1- Column may contain string values and non-float values
// 2- Number of unique values in the column is only a small fraction of
// all values in that column (Lower than or equal to %20 of all values)
$numericValues = array_filter($columnValues, 'is_numeric');
$floatValues = array_filter($columnValues, 'is_float');
if ($floatValues) {
return false;
}
if (count($numericValues) !== $count) {
return true;
}
$distinctValues = array_count_values($columnValues);
return count($distinctValues) <= $count / 5;
}
/**
* This method is used to set number of columns to be used
* when deciding a split at an internal node of the tree. <br>
* If the value is given 0, then all features are used (default behaviour),
* otherwise the given value will be used as a maximum for number of columns
* randomly selected for each split operation.
*
* @param int $numFeatures
* @return $this
* @throws InvalidArgumentException
*/
public function setNumFeatures(int $numFeatures)
{
if ($numFeatures < 0) {
throw new InvalidArgumentException('Selected column count should be greater or equal to zero');
}
$this->numUsableFeatures = $numFeatures;
return $this;
}
/**
* Used to set predefined features to consider while deciding which column to use for a split
*
* @param array $selectedFeatures
*/
protected function setSelectedFeatures(array $selectedFeatures)
{
$this->selectedFeatures = $selectedFeatures;
}
/**
* A string array to represent columns. Useful when HTML output or
* column importances are desired to be inspected.
*
* @param array $names
* @return $this
* @throws InvalidArgumentException
*/
public function setColumnNames(array $names)
{
if ($this->featureCount !== 0 && count($names) !== $this->featureCount) {
throw new InvalidArgumentException(sprintf('Length of the given array should be equal to feature count %s', $this->featureCount));
}
$this->columnNames = $names;
return $this;
}
/**
* @return string
*/
public function getHtml()
{
return $this->tree->getHTML($this->columnNames);
}
/**
* This will return an array including an importance value for
* each column in the given dataset. The importance values are
* normalized and their total makes 1.<br/>
*
* @return array
*/
public function getFeatureImportances()
{
if ($this->featureImportances !== null) {
return $this->featureImportances;
}
$sampleCount = count($this->samples);
$this->featureImportances = [];
foreach ($this->columnNames as $column => $columnName) {
$nodes = $this->getSplitNodesByColumn($column, $this->tree);
$importance = 0;
foreach ($nodes as $node) {
$importance += $node->getNodeImpurityDecrease($sampleCount);
}
$this->featureImportances[$columnName] = $importance;
}
// Normalize & sort the importances
$total = array_sum($this->featureImportances);
if ($total > 0) {
foreach ($this->featureImportances as &$importance) {
$importance /= $total;
}
arsort($this->featureImportances);
}
return $this->featureImportances;
}
/**
* Collects and returns an array of internal nodes that use the given
* column as a split criterion
*
* @param int $column
* @param DecisionTreeLeaf $node
* @return array
*/
protected function getSplitNodesByColumn(int $column, DecisionTreeLeaf $node) : array
{
if (!$node || $node->isTerminal) {
return [];
}
$nodes = [];
if ($node->columnIndex === $column) {
$nodes[] = $node;
}
$lNodes = [];
$rNodes = [];
if ($node->leftLeaf) {
$lNodes = $this->getSplitNodesByColumn($column, $node->leftLeaf);
}
if ($node->rightLeaf) {
$rNodes = $this->getSplitNodesByColumn($column, $node->rightLeaf);
}
$nodes = array_merge($nodes, $lNodes, $rNodes);
return $nodes;
}
/**
* @param array $sample
* @return mixed
*/
protected function predictSample(array $sample)
{
$node = $this->tree;
do {
if ($node->isTerminal) {
break;
}
if ($node->evaluate($sample)) {
$node = $node->leftLeaf;
} else {
$node = $node->rightLeaf;
}
} while ($node);
return $node ? $node->classValue : $this->labels[0];
}
}
@@ -0,0 +1,171 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification\DecisionTree;
class DecisionTreeLeaf
{
/**
* @var string
*/
public $value;
/**
* @var float
*/
public $numericValue;
/**
* @var string
*/
public $operator;
/**
* @var int
*/
public $columnIndex;
/**
* @var DecisionTreeLeaf
*/
public $leftLeaf = null;
/**
* @var DecisionTreeLeaf
*/
public $rightLeaf= null;
/**
* @var array
*/
public $records = [];
/**
* Class value represented by the leaf, this value is non-empty
* only for terminal leaves
*
* @var string
*/
public $classValue = '';
/**
* @var bool
*/
public $isTerminal = false;
/**
* @var bool
*/
public $isContinuous = false;
/**
* @var float
*/
public $giniIndex = 0;
/**
* @var int
*/
public $level = 0;
/**
* @param array $record
* @return bool
*/
public function evaluate($record)
{
$recordField = $record[$this->columnIndex];
if ($this->isContinuous) {
$op = $this->operator;
$value= $this->numericValue;
$recordField = strval($recordField);
eval("\$result = $recordField $op $value;");
return $result;
}
return $recordField == $this->value;
}
/**
* Returns Mean Decrease Impurity (MDI) in the node.
* For terminal nodes, this value is equal to 0
*
* @return float
*/
public function getNodeImpurityDecrease(int $parentRecordCount)
{
if ($this->isTerminal) {
return 0.0;
}
$nodeSampleCount = (float)count($this->records);
$iT = $this->giniIndex;
if ($this->leftLeaf) {
$pL = count($this->leftLeaf->records)/$nodeSampleCount;
$iT -= $pL * $this->leftLeaf->giniIndex;
}
if ($this->rightLeaf) {
$pR = count($this->rightLeaf->records)/$nodeSampleCount;
$iT -= $pR * $this->rightLeaf->giniIndex;
}
return $iT * $nodeSampleCount / $parentRecordCount;
}
/**
* Returns HTML representation of the node including children nodes
*
* @param $columnNames
* @return string
*/
public function getHTML($columnNames = null)
{
if ($this->isTerminal) {
$value = "<b>$this->classValue</b>";
} else {
$value = $this->value;
if ($columnNames !== null) {
$col = $columnNames[$this->columnIndex];
} else {
$col = "col_$this->columnIndex";
}
if (! preg_match("/^[<>=]{1,2}/", $value)) {
$value = "=$value";
}
$value = "<b>$col $value</b><br>Gini: ". number_format($this->giniIndex, 2);
}
$str = "<table ><tr><td colspan=3 align=center style='border:1px solid;'>
$value</td></tr>";
if ($this->leftLeaf || $this->rightLeaf) {
$str .='<tr>';
if ($this->leftLeaf) {
$str .="<td valign=top><b>| Yes</b><br>" . $this->leftLeaf->getHTML($columnNames) . "</td>";
} else {
$str .='<td></td>';
}
$str .='<td>&nbsp;</td>';
if ($this->rightLeaf) {
$str .="<td valign=top align=right><b>No |</b><br>" . $this->rightLeaf->getHTML($columnNames) . "</td>";
} else {
$str .='<td></td>';
}
$str .= '</tr>';
}
$str .= '</table>';
return $str;
}
/**
* HTML representation of the tree without column names
*
* @return string
*/
public function __toString()
{
return $this->getHTML();
}
}
@@ -0,0 +1,265 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification\Ensemble;
use Phpml\Classification\Linear\DecisionStump;
use Phpml\Classification\WeightedClassifier;
use Phpml\Math\Statistic\Mean;
use Phpml\Math\Statistic\StandardDeviation;
use Phpml\Classification\Classifier;
use Phpml\Helper\Predictable;
use Phpml\Helper\Trainable;
class AdaBoost implements Classifier
{
use Predictable, Trainable;
/**
* Actual labels given in the targets array
* @var array
*/
protected $labels = [];
/**
* @var int
*/
protected $sampleCount;
/**
* @var int
*/
protected $featureCount;
/**
* Number of maximum iterations to be done
*
* @var int
*/
protected $maxIterations;
/**
* Sample weights
*
* @var array
*/
protected $weights = [];
/**
* List of selected 'weak' classifiers
*
* @var array
*/
protected $classifiers = [];
/**
* Base classifier weights
*
* @var array
*/
protected $alpha = [];
/**
* @var string
*/
protected $baseClassifier = DecisionStump::class;
/**
* @var array
*/
protected $classifierOptions = [];
/**
* ADAptive BOOSTing (AdaBoost) is an ensemble algorithm to
* improve classification performance of 'weak' classifiers such as
* DecisionStump (default base classifier of AdaBoost).
*
*/
public function __construct(int $maxIterations = 50)
{
$this->maxIterations = $maxIterations;
}
/**
* Sets the base classifier that will be used for boosting (default = DecisionStump)
*
* @param string $baseClassifier
* @param array $classifierOptions
*/
public function setBaseClassifier(string $baseClassifier = DecisionStump::class, array $classifierOptions = [])
{
$this->baseClassifier = $baseClassifier;
$this->classifierOptions = $classifierOptions;
}
/**
* @param array $samples
* @param array $targets
*/
public function train(array $samples, array $targets)
{
// Initialize usual variables
$this->labels = array_keys(array_count_values($targets));
if (count($this->labels) != 2) {
throw new \Exception("AdaBoost is a binary classifier and can classify between two classes only");
}
// Set all target values to either -1 or 1
$this->labels = [1 => $this->labels[0], -1 => $this->labels[1]];
foreach ($targets as $target) {
$this->targets[] = $target == $this->labels[1] ? 1 : -1;
}
$this->samples = array_merge($this->samples, $samples);
$this->featureCount = count($samples[0]);
$this->sampleCount = count($this->samples);
// Initialize AdaBoost parameters
$this->weights = array_fill(0, $this->sampleCount, 1.0 / $this->sampleCount);
$this->classifiers = [];
$this->alpha = [];
// Execute the algorithm for a maximum number of iterations
$currIter = 0;
while ($this->maxIterations > $currIter++) {
// Determine the best 'weak' classifier based on current weights
$classifier = $this->getBestClassifier();
$errorRate = $this->evaluateClassifier($classifier);
// Update alpha & weight values at each iteration
$alpha = $this->calculateAlpha($errorRate);
$this->updateWeights($classifier, $alpha);
$this->classifiers[] = $classifier;
$this->alpha[] = $alpha;
}
}
/**
* Returns the classifier with the lowest error rate with the
* consideration of current sample weights
*
* @return Classifier
*/
protected function getBestClassifier()
{
$ref = new \ReflectionClass($this->baseClassifier);
if ($this->classifierOptions) {
$classifier = $ref->newInstanceArgs($this->classifierOptions);
} else {
$classifier = $ref->newInstance();
}
if (is_subclass_of($classifier, WeightedClassifier::class)) {
$classifier->setSampleWeights($this->weights);
$classifier->train($this->samples, $this->targets);
} else {
list($samples, $targets) = $this->resample();
$classifier->train($samples, $targets);
}
return $classifier;
}
/**
* Resamples the dataset in accordance with the weights and
* returns the new dataset
*
* @return array
*/
protected function resample()
{
$weights = $this->weights;
$std = StandardDeviation::population($weights);
$mean= Mean::arithmetic($weights);
$min = min($weights);
$minZ= (int)round(($min - $mean) / $std);
$samples = [];
$targets = [];
foreach ($weights as $index => $weight) {
$z = (int)round(($weight - $mean) / $std) - $minZ + 1;
for ($i=0; $i < $z; $i++) {
if (rand(0, 1) == 0) {
continue;
}
$samples[] = $this->samples[$index];
$targets[] = $this->targets[$index];
}
}
return [$samples, $targets];
}
/**
* Evaluates the classifier and returns the classification error rate
*
* @param Classifier $classifier
*/
protected function evaluateClassifier(Classifier $classifier)
{
$total = (float) array_sum($this->weights);
$wrong = 0;
foreach ($this->samples as $index => $sample) {
$predicted = $classifier->predict($sample);
if ($predicted != $this->targets[$index]) {
$wrong += $this->weights[$index];
}
}
return $wrong / $total;
}
/**
* Calculates alpha of a classifier
*
* @param float $errorRate
* @return float
*/
protected function calculateAlpha(float $errorRate)
{
if ($errorRate == 0) {
$errorRate = 1e-10;
}
return 0.5 * log((1 - $errorRate) / $errorRate);
}
/**
* Updates the sample weights
*
* @param Classifier $classifier
* @param float $alpha
*/
protected function updateWeights(Classifier $classifier, float $alpha)
{
$sumOfWeights = array_sum($this->weights);
$weightsT1 = [];
foreach ($this->weights as $index => $weight) {
$desired = $this->targets[$index];
$output = $classifier->predict($this->samples[$index]);
$weight *= exp(-$alpha * $desired * $output) / $sumOfWeights;
$weightsT1[] = $weight;
}
$this->weights = $weightsT1;
}
/**
* @param array $sample
* @return mixed
*/
public function predictSample(array $sample)
{
$sum = 0;
foreach ($this->alpha as $index => $alpha) {
$h = $this->classifiers[$index]->predict($sample);
$sum += $h * $alpha;
}
return $this->labels[ $sum > 0 ? 1 : -1];
}
}
@@ -0,0 +1,194 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification\Ensemble;
use Phpml\Helper\Predictable;
use Phpml\Helper\Trainable;
use Phpml\Classification\Classifier;
use Phpml\Classification\DecisionTree;
class Bagging implements Classifier
{
use Trainable, Predictable;
/**
* @var int
*/
protected $numSamples;
/**
* @var array
*/
private $targets = [];
/**
* @var int
*/
protected $featureCount = 0;
/**
* @var int
*/
protected $numClassifier;
/**
* @var Classifier
*/
protected $classifier = DecisionTree::class;
/**
* @var array
*/
protected $classifierOptions = ['depth' => 20];
/**
* @var array
*/
protected $classifiers;
/**
* @var float
*/
protected $subsetRatio = 0.7;
/**
* @var array
*/
private $samples = [];
/**
* Creates an ensemble classifier with given number of base classifiers<br>
* Default number of base classifiers is 100.
* The more number of base classifiers, the better performance but at the cost of procesing time
*
* @param int $numClassifier
*/
public function __construct($numClassifier = 50)
{
$this->numClassifier = $numClassifier;
}
/**
* This method determines the ratio of samples used to create the 'bootstrap' subset,
* e.g., random samples drawn from the original dataset with replacement (allow repeats),
* to train each base classifier.
*
* @param float $ratio
* @return $this
* @throws Exception
*/
public function setSubsetRatio(float $ratio)
{
if ($ratio < 0.1 || $ratio > 1.0) {
throw new \Exception("Subset ratio should be between 0.1 and 1.0");
}
$this->subsetRatio = $ratio;
return $this;
}
/**
* This method is used to set the base classifier. Default value is
* DecisionTree::class, but any class that implements the <i>Classifier</i>
* can be used. <br>
* While giving the parameters of the classifier, the values should be
* given in the order they are in the constructor of the classifier and parameter
* names are neglected.
*
* @param string $classifier
* @param array $classifierOptions
* @return $this
*/
public function setClassifer(string $classifier, array $classifierOptions = [])
{
$this->classifier = $classifier;
$this->classifierOptions = $classifierOptions;
return $this;
}
/**
* @param array $samples
* @param array $targets
*/
public function train(array $samples, array $targets)
{
$this->samples = array_merge($this->samples, $samples);
$this->targets = array_merge($this->targets, $targets);
$this->featureCount = count($samples[0]);
$this->numSamples = count($this->samples);
// Init classifiers and train them with bootstrap samples
$this->classifiers = $this->initClassifiers();
$index = 0;
foreach ($this->classifiers as $classifier) {
list($samples, $targets) = $this->getRandomSubset($index);
$classifier->train($samples, $targets);
++$index;
}
}
/**
* @param int $index
* @return array
*/
protected function getRandomSubset(int $index)
{
$samples = [];
$targets = [];
srand($index);
$bootstrapSize = $this->subsetRatio * $this->numSamples;
for ($i=0; $i < $bootstrapSize; $i++) {
$rand = rand(0, $this->numSamples - 1);
$samples[] = $this->samples[$rand];
$targets[] = $this->targets[$rand];
}
return [$samples, $targets];
}
/**
* @return array
*/
protected function initClassifiers()
{
$classifiers = [];
for ($i=0; $i<$this->numClassifier; $i++) {
$ref = new \ReflectionClass($this->classifier);
if ($this->classifierOptions) {
$obj = $ref->newInstanceArgs($this->classifierOptions);
} else {
$obj = $ref->newInstance();
}
$classifiers[] = $this->initSingleClassifier($obj, $i);
}
return $classifiers;
}
/**
* @param Classifier $classifier
* @param int $index
* @return Classifier
*/
protected function initSingleClassifier($classifier, $index)
{
return $classifier;
}
/**
* @param array $sample
* @return mixed
*/
protected function predictSample(array $sample)
{
$predictions = [];
foreach ($this->classifiers as $classifier) {
/* @var $classifier Classifier */
$predictions[] = $classifier->predict($sample);
}
$counts = array_count_values($predictions);
arsort($counts);
reset($counts);
return key($counts);
}
}
@@ -0,0 +1,153 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification\Ensemble;
use Phpml\Classification\DecisionTree;
use Phpml\Classification\Classifier;
class RandomForest extends Bagging
{
/**
* @var float|string
*/
protected $featureSubsetRatio = 'log';
/**
* @var array
*/
protected $columnNames = null;
/**
* Initializes RandomForest with the given number of trees. More trees
* may increase the prediction performance while it will also substantially
* increase the processing time and the required memory
*
* @param type $numClassifier
*/
public function __construct($numClassifier = 50)
{
parent::__construct($numClassifier);
$this->setSubsetRatio(1.0);
}
/**
* This method is used to determine how many of the original columns (features)
* will be used to construct subsets to train base classifiers.<br>
*
* Allowed values: 'sqrt', 'log' or any float number between 0.1 and 1.0 <br>
*
* Default value for the ratio is 'log' which results in log(numFeatures, 2) + 1
* features to be taken into consideration while selecting subspace of features
*
* @param mixed $ratio string or float should be given
* @return $this
* @throws Exception
*/
public function setFeatureSubsetRatio($ratio)
{
if (is_float($ratio) && ($ratio < 0.1 || $ratio > 1.0)) {
throw new \Exception("When a float given, feature subset ratio should be between 0.1 and 1.0");
}
if (is_string($ratio) && $ratio != 'sqrt' && $ratio != 'log') {
throw new \Exception("When a string given, feature subset ratio can only be 'sqrt' or 'log' ");
}
$this->featureSubsetRatio = $ratio;
return $this;
}
/**
* RandomForest algorithm is usable *only* with DecisionTree
*
* @param string $classifier
* @param array $classifierOptions
* @return $this
*/
public function setClassifer(string $classifier, array $classifierOptions = [])
{
if ($classifier != DecisionTree::class) {
throw new \Exception("RandomForest can only use DecisionTree as base classifier");
}
return parent::setClassifer($classifier, $classifierOptions);
}
/**
* This will return an array including an importance value for
* each column in the given dataset. Importance values for a column
* is the average importance of that column in all trees in the forest
*
* @return array
*/
public function getFeatureImportances()
{
// Traverse each tree and sum importance of the columns
$sum = [];
foreach ($this->classifiers as $tree) {
/* @var $tree DecisionTree */
$importances = $tree->getFeatureImportances();
foreach ($importances as $column => $importance) {
if (array_key_exists($column, $sum)) {
$sum[$column] += $importance;
} else {
$sum[$column] = $importance;
}
}
}
// Normalize & sort the importance values
$total = array_sum($sum);
foreach ($sum as &$importance) {
$importance /= $total;
}
arsort($sum);
return $sum;
}
/**
* A string array to represent the columns is given. They are useful
* when trying to print some information about the trees such as feature importances
*
* @param array $names
* @return $this
*/
public function setColumnNames(array $names)
{
$this->columnNames = $names;
return $this;
}
/**
* @param DecisionTree $classifier
* @param int $index
* @return DecisionTree
*/
protected function initSingleClassifier($classifier, $index)
{
if (is_float($this->featureSubsetRatio)) {
$featureCount = (int)($this->featureSubsetRatio * $this->featureCount);
} elseif ($this->featureCount == 'sqrt') {
$featureCount = (int)sqrt($this->featureCount) + 1;
} else {
$featureCount = (int)log($this->featureCount, 2) + 1;
}
if ($featureCount >= $this->featureCount) {
$featureCount = $this->featureCount;
}
if ($this->columnNames === null) {
$this->columnNames = range(0, $this->featureCount - 1);
}
return $classifier
->setColumnNames($this->columnNames)
->setNumFeatures($featureCount);
}
}
@@ -0,0 +1,82 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification;
use Phpml\Helper\Predictable;
use Phpml\Helper\Trainable;
use Phpml\Math\Distance;
use Phpml\Math\Distance\Euclidean;
class KNearestNeighbors implements Classifier
{
use Trainable, Predictable;
/**
* @var int
*/
private $k;
/**
* @var Distance
*/
private $distanceMetric;
/**
* @param int $k
* @param Distance|null $distanceMetric (if null then Euclidean distance as default)
*/
public function __construct(int $k = 3, Distance $distanceMetric = null)
{
if (null === $distanceMetric) {
$distanceMetric = new Euclidean();
}
$this->k = $k;
$this->samples = [];
$this->targets = [];
$this->distanceMetric = $distanceMetric;
}
/**
* @param array $sample
*
* @return mixed
*/
protected function predictSample(array $sample)
{
$distances = $this->kNeighborsDistances($sample);
$predictions = array_combine(array_values($this->targets), array_fill(0, count($this->targets), 0));
foreach ($distances as $index => $distance) {
++$predictions[$this->targets[$index]];
}
arsort($predictions);
reset($predictions);
return key($predictions);
}
/**
* @param array $sample
*
* @return array
*
* @throws \Phpml\Exception\InvalidArgumentException
*/
private function kNeighborsDistances(array $sample)
{
$distances = [];
foreach ($this->samples as $index => $neighbor) {
$distances[$index] = $this->distanceMetric->distance($sample, $neighbor);
}
asort($distances);
return array_slice($distances, 0, $this->k, true);
}
}
@@ -0,0 +1,77 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification\Linear;
use Phpml\Classification\Classifier;
class Adaline extends Perceptron
{
/**
* Batch training is the default Adaline training algorithm
*/
const BATCH_TRAINING = 1;
/**
* Online training: Stochastic gradient descent learning
*/
const ONLINE_TRAINING = 2;
/**
* Training type may be either 'Batch' or 'Online' learning
*
* @var string
*/
protected $trainingType;
/**
* Initalize an Adaline (ADAptive LInear NEuron) classifier with given learning rate and maximum
* number of iterations used while training the classifier <br>
*
* Learning rate should be a float value between 0.0(exclusive) and 1.0 (inclusive) <br>
* Maximum number of iterations can be an integer value greater than 0 <br>
* If normalizeInputs is set to true, then every input given to the algorithm will be standardized
* by use of standard deviation and mean calculation
*
* @param int $learningRate
* @param int $maxIterations
*/
public function __construct(float $learningRate = 0.001, int $maxIterations = 1000,
bool $normalizeInputs = true, int $trainingType = self::BATCH_TRAINING)
{
if (! in_array($trainingType, [self::BATCH_TRAINING, self::ONLINE_TRAINING])) {
throw new \Exception("Adaline can only be trained with batch and online/stochastic gradient descent algorithm");
}
$this->trainingType = $trainingType;
parent::__construct($learningRate, $maxIterations, $normalizeInputs);
}
/**
* Adapts the weights with respect to given samples and targets
* by use of gradient descent learning rule
*
* @param array $samples
* @param array $targets
*/
protected function runTraining(array $samples, array $targets)
{
// The cost function is the sum of squares
$callback = function ($weights, $sample, $target) {
$this->weights = $weights;
$output = $this->output($sample);
$gradient = $output - $target;
$error = $gradient ** 2;
return [$error, $gradient];
};
$isBatch = $this->trainingType == self::BATCH_TRAINING;
return parent::runGradientDescent($samples, $targets, $callback, $isBatch);
}
}
@@ -0,0 +1,362 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification\Linear;
use Phpml\Helper\Predictable;
use Phpml\Helper\OneVsRest;
use Phpml\Classification\WeightedClassifier;
use Phpml\Classification\DecisionTree;
class DecisionStump extends WeightedClassifier
{
use Predictable, OneVsRest;
const AUTO_SELECT = -1;
/**
* @var int
*/
protected $givenColumnIndex;
/**
* @var array
*/
protected $binaryLabels;
/**
* Lowest error rate obtained while training/optimizing the model
*
* @var float
*/
protected $trainingErrorRate;
/**
* @var int
*/
protected $column;
/**
* @var mixed
*/
protected $value;
/**
* @var string
*/
protected $operator;
/**
* @var array
*/
protected $columnTypes;
/**
* @var int
*/
protected $featureCount;
/**
* @var float
*/
protected $numSplitCount = 100.0;
/**
* Distribution of samples in the leaves
*
* @var array
*/
protected $prob;
/**
* A DecisionStump classifier is a one-level deep DecisionTree. It is generally
* used with ensemble algorithms as in the weak classifier role. <br>
*
* If columnIndex is given, then the stump tries to produce a decision node
* on this column, otherwise in cases given the value of -1, the stump itself
* decides which column to take for the decision (Default DecisionTree behaviour)
*
* @param int $columnIndex
*/
public function __construct(int $columnIndex = self::AUTO_SELECT)
{
$this->givenColumnIndex = $columnIndex;
}
/**
* @param array $samples
* @param array $targets
* @throws \Exception
*/
protected function trainBinary(array $samples, array $targets, array $labels)
{
$this->binaryLabels = $labels;
$this->featureCount = count($samples[0]);
// If a column index is given, it should be among the existing columns
if ($this->givenColumnIndex > count($samples[0]) - 1) {
$this->givenColumnIndex = self::AUTO_SELECT;
}
// Check the size of the weights given.
// If none given, then assign 1 as a weight to each sample
if ($this->weights) {
$numWeights = count($this->weights);
if ($numWeights != count($samples)) {
throw new \Exception("Number of sample weights does not match with number of samples");
}
} else {
$this->weights = array_fill(0, count($samples), 1);
}
// Determine type of each column as either "continuous" or "nominal"
$this->columnTypes = DecisionTree::getColumnTypes($samples);
// Try to find the best split in the columns of the dataset
// by calculating error rate for each split point in each column
$columns = range(0, count($samples[0]) - 1);
if ($this->givenColumnIndex != self::AUTO_SELECT) {
$columns = [$this->givenColumnIndex];
}
$bestSplit = [
'value' => 0, 'operator' => '',
'prob' => [], 'column' => 0,
'trainingErrorRate' => 1.0];
foreach ($columns as $col) {
if ($this->columnTypes[$col] == DecisionTree::CONTINUOUS) {
$split = $this->getBestNumericalSplit($samples, $targets, $col);
} else {
$split = $this->getBestNominalSplit($samples, $targets, $col);
}
if ($split['trainingErrorRate'] < $bestSplit['trainingErrorRate']) {
$bestSplit = $split;
}
}
// Assign determined best values to the stump
foreach ($bestSplit as $name => $value) {
$this->{$name} = $value;
}
}
/**
* While finding best split point for a numerical valued column,
* DecisionStump looks for equally distanced values between minimum and maximum
* values in the column. Given <i>$count</i> value determines how many split
* points to be probed. The more split counts, the better performance but
* worse processing time (Default value is 10.0)
*
* @param float $count
*/
public function setNumericalSplitCount(float $count)
{
$this->numSplitCount = $count;
}
/**
* Determines best split point for the given column
*
* @param array $samples
* @param array $targets
* @param int $col
*
* @return array
*/
protected function getBestNumericalSplit(array $samples, array $targets, int $col)
{
$values = array_column($samples, $col);
// Trying all possible points may be accomplished in two general ways:
// 1- Try all values in the $samples array ($values)
// 2- Artificially split the range of values into several parts and try them
// We choose the second one because it is faster in larger datasets
$minValue = min($values);
$maxValue = max($values);
$stepSize = ($maxValue - $minValue) / $this->numSplitCount;
$split = null;
foreach (['<=', '>'] as $operator) {
// Before trying all possible split points, let's first try
// the average value for the cut point
$threshold = array_sum($values) / (float) count($values);
list($errorRate, $prob) = $this->calculateErrorRate($targets, $threshold, $operator, $values);
if ($split == null || $errorRate < $split['trainingErrorRate']) {
$split = ['value' => $threshold, 'operator' => $operator,
'prob' => $prob, 'column' => $col,
'trainingErrorRate' => $errorRate];
}
// Try other possible points one by one
for ($step = $minValue; $step <= $maxValue; $step+= $stepSize) {
$threshold = (float)$step;
list($errorRate, $prob) = $this->calculateErrorRate($targets, $threshold, $operator, $values);
if ($errorRate < $split['trainingErrorRate']) {
$split = ['value' => $threshold, 'operator' => $operator,
'prob' => $prob, 'column' => $col,
'trainingErrorRate' => $errorRate];
}
}// for
}
return $split;
}
/**
* @param array $samples
* @param array $targets
* @param int $col
*
* @return array
*/
protected function getBestNominalSplit(array $samples, array $targets, int $col) : array
{
$values = array_column($samples, $col);
$valueCounts = array_count_values($values);
$distinctVals= array_keys($valueCounts);
$split = null;
foreach (['=', '!='] as $operator) {
foreach ($distinctVals as $val) {
list($errorRate, $prob) = $this->calculateErrorRate($targets, $val, $operator, $values);
if ($split == null || $split['trainingErrorRate'] < $errorRate) {
$split = ['value' => $val, 'operator' => $operator,
'prob' => $prob, 'column' => $col,
'trainingErrorRate' => $errorRate];
}
}
}
return $split;
}
/**
*
* @param type $leftValue
* @param type $operator
* @param type $rightValue
*
* @return boolean
*/
protected function evaluate($leftValue, $operator, $rightValue)
{
switch ($operator) {
case '>': return $leftValue > $rightValue;
case '>=': return $leftValue >= $rightValue;
case '<': return $leftValue < $rightValue;
case '<=': return $leftValue <= $rightValue;
case '=': return $leftValue === $rightValue;
case '!=':
case '<>': return $leftValue !== $rightValue;
}
return false;
}
/**
* Calculates the ratio of wrong predictions based on the new threshold
* value given as the parameter
*
* @param array $targets
* @param float $threshold
* @param string $operator
* @param array $values
*
* @return array
*/
protected function calculateErrorRate(array $targets, float $threshold, string $operator, array $values) : array
{
$wrong = 0.0;
$prob = [];
$leftLabel = $this->binaryLabels[0];
$rightLabel= $this->binaryLabels[1];
foreach ($values as $index => $value) {
if ($this->evaluate($value, $operator, $threshold)) {
$predicted = $leftLabel;
} else {
$predicted = $rightLabel;
}
$target = $targets[$index];
if (strval($predicted) != strval($targets[$index])) {
$wrong += $this->weights[$index];
}
if (! isset($prob[$predicted][$target])) {
$prob[$predicted][$target] = 0;
}
$prob[$predicted][$target]++;
}
// Calculate probabilities: Proportion of labels in each leaf
$dist = array_combine($this->binaryLabels, array_fill(0, 2, 0.0));
foreach ($prob as $leaf => $counts) {
$leafTotal = (float)array_sum($prob[$leaf]);
foreach ($counts as $label => $count) {
if (strval($leaf) == strval($label)) {
$dist[$leaf] = $count / $leafTotal;
}
}
}
return [$wrong / (float) array_sum($this->weights), $dist];
}
/**
* Returns the probability of the sample of belonging to the given label
*
* Probability of a sample is calculated as the proportion of the label
* within the labels of the training samples in the decision node
*
* @param array $sample
* @param mixed $label
*
* @return float
*/
protected function predictProbability(array $sample, $label) : float
{
$predicted = $this->predictSampleBinary($sample);
if (strval($predicted) == strval($label)) {
return $this->prob[$label];
}
return 0.0;
}
/**
* @param array $sample
*
* @return mixed
*/
protected function predictSampleBinary(array $sample)
{
if ($this->evaluate($sample[$this->column], $this->operator, $this->value)) {
return $this->binaryLabels[0];
}
return $this->binaryLabels[1];
}
/**
* @return void
*/
protected function resetBinary()
{
}
/**
* @return string
*/
public function __toString()
{
return "IF $this->column $this->operator $this->value " .
"THEN " . $this->binaryLabels[0] . " ".
"ELSE " . $this->binaryLabels[1];
}
}
@@ -0,0 +1,281 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification\Linear;
use Phpml\Classification\Classifier;
use Phpml\Helper\Optimizer\ConjugateGradient;
class LogisticRegression extends Adaline
{
/**
* Batch training: Gradient descent algorithm (default)
*/
const BATCH_TRAINING = 1;
/**
* Online training: Stochastic gradient descent learning
*/
const ONLINE_TRAINING = 2;
/**
* Conjugate Batch: Conjugate Gradient algorithm
*/
const CONJUGATE_GRAD_TRAINING = 3;
/**
* Cost function to optimize: 'log' and 'sse' are supported <br>
* - 'log' : log likelihood <br>
* - 'sse' : sum of squared errors <br>
*
* @var string
*/
protected $costFunction = 'sse';
/**
* Regularization term: only 'L2' is supported
*
* @var string
*/
protected $penalty = 'L2';
/**
* Lambda (λ) parameter of regularization term. If λ is set to 0, then
* regularization term is cancelled.
*
* @var float
*/
protected $lambda = 0.5;
/**
* Initalize a Logistic Regression classifier with maximum number of iterations
* and learning rule to be applied <br>
*
* Maximum number of iterations can be an integer value greater than 0 <br>
* If normalizeInputs is set to true, then every input given to the algorithm will be standardized
* by use of standard deviation and mean calculation <br>
*
* Cost function can be 'log' for log-likelihood and 'sse' for sum of squared errors <br>
*
* Penalty (Regularization term) can be 'L2' or empty string to cancel penalty term
*
* @param int $maxIterations
* @param bool $normalizeInputs
* @param int $trainingType
* @param string $cost
* @param string $penalty
*
* @throws \Exception
*/
public function __construct(int $maxIterations = 500, bool $normalizeInputs = true,
int $trainingType = self::CONJUGATE_GRAD_TRAINING, string $cost = 'sse',
string $penalty = 'L2')
{
$trainingTypes = range(self::BATCH_TRAINING, self::CONJUGATE_GRAD_TRAINING);
if (! in_array($trainingType, $trainingTypes)) {
throw new \Exception("Logistic regression can only be trained with " .
"batch (gradient descent), online (stochastic gradient descent) " .
"or conjugate batch (conjugate gradients) algorithms");
}
if (! in_array($cost, ['log', 'sse'])) {
throw new \Exception("Logistic regression cost function can be one of the following: \n" .
"'log' for log-likelihood and 'sse' for sum of squared errors");
}
if ($penalty != '' && strtoupper($penalty) !== 'L2') {
throw new \Exception("Logistic regression supports only 'L2' regularization");
}
$this->learningRate = 0.001;
parent::__construct($this->learningRate, $maxIterations, $normalizeInputs);
$this->trainingType = $trainingType;
$this->costFunction = $cost;
$this->penalty = $penalty;
}
/**
* Sets the learning rate if gradient descent algorithm is
* selected for training
*
* @param float $learningRate
*/
public function setLearningRate(float $learningRate)
{
$this->learningRate = $learningRate;
}
/**
* Lambda (λ) parameter of regularization term. If 0 is given,
* then the regularization term is cancelled
*
* @param float $lambda
*/
public function setLambda(float $lambda)
{
$this->lambda = $lambda;
}
/**
* Adapts the weights with respect to given samples and targets
* by use of selected solver
*
* @param array $samples
* @param array $targets
*/
protected function runTraining(array $samples, array $targets)
{
$callback = $this->getCostFunction();
switch ($this->trainingType) {
case self::BATCH_TRAINING:
return $this->runGradientDescent($samples, $targets, $callback, true);
case self::ONLINE_TRAINING:
return $this->runGradientDescent($samples, $targets, $callback, false);
case self::CONJUGATE_GRAD_TRAINING:
return $this->runConjugateGradient($samples, $targets, $callback);
}
}
/**
* Executes Conjugate Gradient method to optimize the
* weights of the LogReg model
*/
protected function runConjugateGradient(array $samples, array $targets, \Closure $gradientFunc)
{
if (empty($this->optimizer)) {
$this->optimizer = (new ConjugateGradient($this->featureCount))
->setMaxIterations($this->maxIterations);
}
$this->weights = $this->optimizer->runOptimization($samples, $targets, $gradientFunc);
$this->costValues = $this->optimizer->getCostValues();
}
/**
* Returns the appropriate callback function for the selected cost function
*
* @return \Closure
*/
protected function getCostFunction()
{
$penalty = 0;
if ($this->penalty == 'L2') {
$penalty = $this->lambda;
}
switch ($this->costFunction) {
case 'log':
/*
* Negative of Log-likelihood cost function to be minimized:
* J(x) = ∑( - y . log(h(x)) - (1 - y) . log(1 - h(x)))
*
* If regularization term is given, then it will be added to the cost:
* for L2 : J(x) = J(x) + λ/m . w
*
* The gradient of the cost function to be used with gradient descent:
* ∇J(x) = -(y - h(x)) = (h(x) - y)
*/
$callback = function ($weights, $sample, $y) use ($penalty) {
$this->weights = $weights;
$hX = $this->output($sample);
// In cases where $hX = 1 or $hX = 0, the log-likelihood
// value will give a NaN, so we fix these values
if ($hX == 1) {
$hX = 1 - 1e-10;
}
if ($hX == 0) {
$hX = 1e-10;
}
$error = -$y * log($hX) - (1 - $y) * log(1 - $hX);
$gradient = $hX - $y;
return [$error, $gradient, $penalty];
};
return $callback;
case 'sse':
/**
* Sum of squared errors or least squared errors cost function:
* J(x) = ∑ (y - h(x))^2
*
* If regularization term is given, then it will be added to the cost:
* for L2 : J(x) = J(x) + λ/m . w
*
* The gradient of the cost function:
* ∇J(x) = -(h(x) - y) . h(x) . (1 - h(x))
*/
$callback = function ($weights, $sample, $y) use ($penalty) {
$this->weights = $weights;
$hX = $this->output($sample);
$error = ($y - $hX) ** 2;
$gradient = -($y - $hX) * $hX * (1 - $hX);
return [$error, $gradient, $penalty];
};
return $callback;
}
}
/**
* Returns the output of the network, a float value between 0.0 and 1.0
*
* @param array $sample
*
* @return float
*/
protected function output(array $sample)
{
$sum = parent::output($sample);
return 1.0 / (1.0 + exp(-$sum));
}
/**
* Returns the class value (either -1 or 1) for the given input
*
* @param array $sample
* @return int
*/
protected function outputClass(array $sample)
{
$output = $this->output($sample);
if (round($output) > 0.5) {
return 1;
}
return -1;
}
/**
* Returns the probability of the sample of belonging to the given label.
*
* The probability is simply taken as the distance of the sample
* to the decision plane.
*
* @param array $sample
* @param mixed $label
*/
protected function predictProbability(array $sample, $label)
{
$predicted = $this->predictSampleBinary($sample);
if (strval($predicted) == strval($label)) {
$sample = $this->checkNormalizedSample($sample);
return abs($this->output($sample) - 0.5);
}
return 0.0;
}
}
@@ -0,0 +1,290 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification\Linear;
use Phpml\Helper\Predictable;
use Phpml\Helper\OneVsRest;
use Phpml\Helper\Optimizer\StochasticGD;
use Phpml\Helper\Optimizer\GD;
use Phpml\Classification\Classifier;
use Phpml\Preprocessing\Normalizer;
use Phpml\IncrementalEstimator;
class Perceptron implements Classifier, IncrementalEstimator
{
use Predictable, OneVsRest;
/**
* @var \Phpml\Helper\Optimizer\Optimizer
*/
protected $optimizer;
/**
* @var array
*/
protected $labels = [];
/**
* @var int
*/
protected $featureCount = 0;
/**
* @var array
*/
protected $weights;
/**
* @var float
*/
protected $learningRate;
/**
* @var int
*/
protected $maxIterations;
/**
* @var Normalizer
*/
protected $normalizer;
/**
* @var bool
*/
protected $enableEarlyStop = true;
/**
* @var array
*/
protected $costValues = [];
/**
* Initalize a perceptron classifier with given learning rate and maximum
* number of iterations used while training the perceptron <br>
*
* Learning rate should be a float value between 0.0(exclusive) and 1.0(inclusive) <br>
* Maximum number of iterations can be an integer value greater than 0
* @param int $learningRate
* @param int $maxIterations
*/
public function __construct(float $learningRate = 0.001, int $maxIterations = 1000,
bool $normalizeInputs = true)
{
if ($learningRate <= 0.0 || $learningRate > 1.0) {
throw new \Exception("Learning rate should be a float value between 0.0(exclusive) and 1.0(inclusive)");
}
if ($maxIterations <= 0) {
throw new \Exception("Maximum number of iterations should be an integer greater than 0");
}
if ($normalizeInputs) {
$this->normalizer = new Normalizer(Normalizer::NORM_STD);
}
$this->learningRate = $learningRate;
$this->maxIterations = $maxIterations;
}
/**
* @param array $samples
* @param array $targets
* @param array $labels
*/
public function partialTrain(array $samples, array $targets, array $labels = [])
{
return $this->trainByLabel($samples, $targets, $labels);
}
/**
* @param array $samples
* @param array $targets
* @param array $labels
*/
public function trainBinary(array $samples, array $targets, array $labels)
{
if ($this->normalizer) {
$this->normalizer->transform($samples);
}
// Set all target values to either -1 or 1
$this->labels = [1 => $labels[0], -1 => $labels[1]];
foreach ($targets as $key => $target) {
$targets[$key] = strval($target) == strval($this->labels[1]) ? 1 : -1;
}
// Set samples and feature count vars
$this->featureCount = count($samples[0]);
$this->runTraining($samples, $targets);
}
protected function resetBinary()
{
$this->labels = [];
$this->optimizer = null;
$this->featureCount = 0;
$this->weights = null;
$this->costValues = [];
}
/**
* Normally enabling early stopping for the optimization procedure may
* help saving processing time while in some cases it may result in
* premature convergence.<br>
*
* If "false" is given, the optimization procedure will always be executed
* for $maxIterations times
*
* @param bool $enable
*/
public function setEarlyStop(bool $enable = true)
{
$this->enableEarlyStop = $enable;
return $this;
}
/**
* Returns the cost values obtained during the training.
*
* @return array
*/
public function getCostValues()
{
return $this->costValues;
}
/**
* Trains the perceptron model with Stochastic Gradient Descent optimization
* to get the correct set of weights
*
* @param array $samples
* @param array $targets
*/
protected function runTraining(array $samples, array $targets)
{
// The cost function is the sum of squares
$callback = function ($weights, $sample, $target) {
$this->weights = $weights;
$prediction = $this->outputClass($sample);
$gradient = $prediction - $target;
$error = $gradient**2;
return [$error, $gradient];
};
$this->runGradientDescent($samples, $targets, $callback);
}
/**
* Executes a Gradient Descent algorithm for
* the given cost function
*
* @param array $samples
* @param array $targets
*/
protected function runGradientDescent(array $samples, array $targets, \Closure $gradientFunc, bool $isBatch = false)
{
$class = $isBatch ? GD::class : StochasticGD::class;
if (empty($this->optimizer)) {
$this->optimizer = (new $class($this->featureCount))
->setLearningRate($this->learningRate)
->setMaxIterations($this->maxIterations)
->setChangeThreshold(1e-6)
->setEarlyStop($this->enableEarlyStop);
}
$this->weights = $this->optimizer->runOptimization($samples, $targets, $gradientFunc);
$this->costValues = $this->optimizer->getCostValues();
}
/**
* Checks if the sample should be normalized and if so, returns the
* normalized sample
*
* @param array $sample
*
* @return array
*/
protected function checkNormalizedSample(array $sample)
{
if ($this->normalizer) {
$samples = [$sample];
$this->normalizer->transform($samples);
$sample = $samples[0];
}
return $sample;
}
/**
* Calculates net output of the network as a float value for the given input
*
* @param array $sample
* @return int
*/
protected function output(array $sample)
{
$sum = 0;
foreach ($this->weights as $index => $w) {
if ($index == 0) {
$sum += $w;
} else {
$sum += $w * $sample[$index - 1];
}
}
return $sum;
}
/**
* Returns the class value (either -1 or 1) for the given input
*
* @param array $sample
* @return int
*/
protected function outputClass(array $sample)
{
return $this->output($sample) > 0 ? 1 : -1;
}
/**
* Returns the probability of the sample of belonging to the given label.
*
* The probability is simply taken as the distance of the sample
* to the decision plane.
*
* @param array $sample
* @param mixed $label
*/
protected function predictProbability(array $sample, $label)
{
$predicted = $this->predictSampleBinary($sample);
if (strval($predicted) == strval($label)) {
$sample = $this->checkNormalizedSample($sample);
return abs($this->output($sample));
}
return 0.0;
}
/**
* @param array $sample
* @return mixed
*/
protected function predictSampleBinary(array $sample)
{
$sample = $this->checkNormalizedSample($sample);
$predictedClass = $this->outputClass($sample);
return $this->labels[ $predictedClass ];
}
}
@@ -0,0 +1,184 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification;
use Phpml\Helper\Predictable;
use Phpml\Helper\Trainable;
use Phpml\Math\Statistic\Mean;
use Phpml\Math\Statistic\StandardDeviation;
class NaiveBayes implements Classifier
{
use Trainable, Predictable;
const CONTINUOS = 1;
const NOMINAL = 2;
const EPSILON = 1e-10;
/**
* @var array
*/
private $std = [];
/**
* @var array
*/
private $mean= [];
/**
* @var array
*/
private $discreteProb = [];
/**
* @var array
*/
private $dataType = [];
/**
* @var array
*/
private $p = [];
/**
* @var int
*/
private $sampleCount = 0;
/**
* @var int
*/
private $featureCount = 0;
/**
* @var array
*/
private $labels = [];
/**
* @param array $samples
* @param array $targets
*/
public function train(array $samples, array $targets)
{
$this->samples = array_merge($this->samples, $samples);
$this->targets = array_merge($this->targets, $targets);
$this->sampleCount = count($this->samples);
$this->featureCount = count($this->samples[0]);
$labelCounts = array_count_values($this->targets);
$this->labels = array_keys($labelCounts);
foreach ($this->labels as $label) {
$samples = $this->getSamplesByLabel($label);
$this->p[$label] = count($samples) / $this->sampleCount;
$this->calculateStatistics($label, $samples);
}
}
/**
* Calculates vital statistics for each label & feature. Stores these
* values in private array in order to avoid repeated calculation
* @param string $label
* @param array $samples
*/
private function calculateStatistics($label, $samples)
{
$this->std[$label] = array_fill(0, $this->featureCount, 0);
$this->mean[$label]= array_fill(0, $this->featureCount, 0);
$this->dataType[$label] = array_fill(0, $this->featureCount, self::CONTINUOS);
$this->discreteProb[$label] = array_fill(0, $this->featureCount, self::CONTINUOS);
for ($i=0; $i<$this->featureCount; $i++) {
// Get the values of nth column in the samples array
// Mean::arithmetic is called twice, can be optimized
$values = array_column($samples, $i);
$numValues = count($values);
// if the values contain non-numeric data,
// then it should be treated as nominal/categorical/discrete column
if ($values !== array_filter($values, 'is_numeric')) {
$this->dataType[$label][$i] = self::NOMINAL;
$this->discreteProb[$label][$i] = array_count_values($values);
$db = &$this->discreteProb[$label][$i];
$db = array_map(function ($el) use ($numValues) {
return $el / $numValues;
}, $db);
} else {
$this->mean[$label][$i] = Mean::arithmetic($values);
// Add epsilon in order to avoid zero stdev
$this->std[$label][$i] = 1e-10 + StandardDeviation::population($values, false);
}
}
}
/**
* Calculates the probability P(label|sample_n)
*
* @param array $sample
* @param int $feature
* @param string $label
* @return float
*/
private function sampleProbability($sample, $feature, $label)
{
$value = $sample[$feature];
if ($this->dataType[$label][$feature] == self::NOMINAL) {
if (! isset($this->discreteProb[$label][$feature][$value]) ||
$this->discreteProb[$label][$feature][$value] == 0) {
return self::EPSILON;
}
return $this->discreteProb[$label][$feature][$value];
}
$std = $this->std[$label][$feature] ;
$mean= $this->mean[$label][$feature];
// Calculate the probability density by use of normal/Gaussian distribution
// Ref: https://en.wikipedia.org/wiki/Normal_distribution
//
// In order to avoid numerical errors because of small or zero values,
// some libraries adopt taking log of calculations such as
// scikit-learn did.
// (See : https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/naive_bayes.py)
$pdf = -0.5 * log(2.0 * pi() * $std * $std);
$pdf -= 0.5 * pow($value - $mean, 2) / ($std * $std);
return $pdf;
}
/**
* Return samples belonging to specific label
* @param string $label
* @return array
*/
private function getSamplesByLabel($label)
{
$samples = [];
for ($i=0; $i<$this->sampleCount; $i++) {
if ($this->targets[$i] == $label) {
$samples[] = $this->samples[$i];
}
}
return $samples;
}
/**
* @param array $sample
* @return mixed
*/
protected function predictSample(array $sample)
{
// Use NaiveBayes assumption for each label using:
// P(label|features) = P(label) * P(feature0|label) * P(feature1|label) .... P(featureN|label)
// Then compare probability for each class to determine which label is most likely
$predictions = [];
foreach ($this->labels as $label) {
$p = $this->p[$label];
for ($i=0; $i<$this->featureCount; $i++) {
$Plf = $this->sampleProbability($sample, $i, $label);
$p += $Plf;
}
$predictions[$label] = $p;
}
arsort($predictions, SORT_NUMERIC);
reset($predictions);
return key($predictions);
}
}
@@ -0,0 +1,31 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification;
use Phpml\SupportVectorMachine\Kernel;
use Phpml\SupportVectorMachine\SupportVectorMachine;
use Phpml\SupportVectorMachine\Type;
class SVC extends SupportVectorMachine implements Classifier
{
/**
* @param int $kernel
* @param float $cost
* @param int $degree
* @param float|null $gamma
* @param float $coef0
* @param float $tolerance
* @param int $cacheSize
* @param bool $shrinking
* @param bool $probabilityEstimates
*/
public function __construct(
int $kernel = Kernel::LINEAR, float $cost = 1.0, int $degree = 3, float $gamma = null, float $coef0 = 0.0,
float $tolerance = 0.001, int $cacheSize = 100, bool $shrinking = true,
bool $probabilityEstimates = false
) {
parent::__construct(Type::C_SVC, $kernel, $cost, 0.5, $degree, $gamma, $coef0, 0.1, $tolerance, $cacheSize, $shrinking, $probabilityEstimates);
}
}
@@ -0,0 +1,23 @@
<?php
declare(strict_types=1);
namespace Phpml\Classification;
abstract class WeightedClassifier implements Classifier
{
/**
* @var array
*/
protected $weights;
/**
* Sets the array including a weight for each sample
*
* @param array $weights
*/
public function setSampleWeights(array $weights)
{
$this->weights = $weights;
}
}
@@ -0,0 +1,15 @@
<?php
declare(strict_types=1);
namespace Phpml\Clustering;
interface Clusterer
{
/**
* @param array $samples
*
* @return array
*/
public function cluster(array $samples);
}
@@ -0,0 +1,111 @@
<?php
declare(strict_types=1);
namespace Phpml\Clustering;
use Phpml\Math\Distance;
use Phpml\Math\Distance\Euclidean;
class DBSCAN implements Clusterer
{
/**
* @var float
*/
private $epsilon;
/**
* @var int
*/
private $minSamples;
/**
* @var Distance
*/
private $distanceMetric;
/**
* @param float $epsilon
* @param int $minSamples
* @param Distance $distanceMetric
*/
public function __construct($epsilon = 0.5, $minSamples = 3, Distance $distanceMetric = null)
{
if (null === $distanceMetric) {
$distanceMetric = new Euclidean();
}
$this->epsilon = $epsilon;
$this->minSamples = $minSamples;
$this->distanceMetric = $distanceMetric;
}
/**
* @param array $samples
*
* @return array
*/
public function cluster(array $samples)
{
$clusters = [];
$visited = [];
foreach ($samples as $index => $sample) {
if (isset($visited[$index])) {
continue;
}
$visited[$index] = true;
$regionSamples = $this->getSamplesInRegion($sample, $samples);
if (count($regionSamples) >= $this->minSamples) {
$clusters[] = $this->expandCluster($regionSamples, $visited);
}
}
return $clusters;
}
/**
* @param array $localSample
* @param array $samples
*
* @return array
*/
private function getSamplesInRegion($localSample, $samples)
{
$region = [];
foreach ($samples as $index => $sample) {
if ($this->distanceMetric->distance($localSample, $sample) < $this->epsilon) {
$region[$index] = $sample;
}
}
return $region;
}
/**
* @param array $samples
* @param array $visited
*
* @return array
*/
private function expandCluster($samples, &$visited)
{
$cluster = [];
foreach ($samples as $index => $sample) {
if (!isset($visited[$index])) {
$visited[$index] = true;
$regionSamples = $this->getSamplesInRegion($sample, $samples);
if (count($regionSamples) > $this->minSamples) {
$cluster = array_merge($regionSamples, $cluster);
}
}
$cluster[] = $sample;
}
return $cluster;
}
}
@@ -0,0 +1,243 @@
<?php
declare(strict_types=1);
namespace Phpml\Clustering;
use Phpml\Clustering\KMeans\Point;
use Phpml\Clustering\KMeans\Cluster;
use Phpml\Clustering\KMeans\Space;
use Phpml\Math\Distance\Euclidean;
class FuzzyCMeans implements Clusterer
{
/**
* @var int
*/
private $clustersNumber;
/**
* @var array|Cluster[]
*/
private $clusters = null;
/**
* @var Space
*/
private $space;
/**
* @var array|float[][]
*/
private $membership;
/**
* @var float
*/
private $fuzziness;
/**
* @var float
*/
private $epsilon;
/**
* @var int
*/
private $maxIterations;
/**
* @var int
*/
private $sampleCount;
/**
* @var array
*/
private $samples;
/**
* @param int $clustersNumber
*
* @throws InvalidArgumentException
*/
public function __construct(int $clustersNumber, float $fuzziness = 2.0, float $epsilon = 1e-2, int $maxIterations = 100)
{
if ($clustersNumber <= 0) {
throw InvalidArgumentException::invalidClustersNumber();
}
$this->clustersNumber = $clustersNumber;
$this->fuzziness = $fuzziness;
$this->epsilon = $epsilon;
$this->maxIterations = $maxIterations;
}
protected function initClusters()
{
// Membership array is a matrix of cluster number by sample counts
// We initilize the membership array with random values
$dim = $this->space->getDimension();
$this->generateRandomMembership($dim, $this->sampleCount);
$this->updateClusters();
}
/**
* @param int $rows
* @param int $cols
*/
protected function generateRandomMembership(int $rows, int $cols)
{
$this->membership = [];
for ($i=0; $i < $rows; $i++) {
$row = [];
$total = 0.0;
for ($k=0; $k < $cols; $k++) {
$val = rand(1, 5) / 10.0;
$row[] = $val;
$total += $val;
}
$this->membership[] = array_map(function ($val) use ($total) {
return $val / $total;
}, $row);
}
}
protected function updateClusters()
{
$dim = $this->space->getDimension();
if (! $this->clusters) {
$this->clusters = [];
for ($i=0; $i<$this->clustersNumber; $i++) {
$this->clusters[] = new Cluster($this->space, array_fill(0, $dim, 0.0));
}
}
for ($i=0; $i<$this->clustersNumber; $i++) {
$cluster = $this->clusters[$i];
$center = $cluster->getCoordinates();
for ($k=0; $k<$dim; $k++) {
$a = $this->getMembershipRowTotal($i, $k, true);
$b = $this->getMembershipRowTotal($i, $k, false);
$center[$k] = $a / $b;
}
$cluster->setCoordinates($center);
}
}
protected function getMembershipRowTotal(int $row, int $col, bool $multiply)
{
$sum = 0.0;
for ($k = 0; $k < $this->sampleCount; $k++) {
$val = pow($this->membership[$row][$k], $this->fuzziness);
if ($multiply) {
$val *= $this->samples[$k][$col];
}
$sum += $val;
}
return $sum;
}
protected function updateMembershipMatrix()
{
for ($i = 0; $i < $this->clustersNumber; $i++) {
for ($k = 0; $k < $this->sampleCount; $k++) {
$distCalc = $this->getDistanceCalc($i, $k);
$this->membership[$i][$k] = 1.0 / $distCalc;
}
}
}
/**
*
* @param int $row
* @param int $col
* @return float
*/
protected function getDistanceCalc(int $row, int $col)
{
$sum = 0.0;
$distance = new Euclidean();
$dist1 = $distance->distance(
$this->clusters[$row]->getCoordinates(),
$this->samples[$col]);
for ($j = 0; $j < $this->clustersNumber; $j++) {
$dist2 = $distance->distance(
$this->clusters[$j]->getCoordinates(),
$this->samples[$col]);
$val = pow($dist1 / $dist2, 2.0 / ($this->fuzziness - 1));
$sum += $val;
}
return $sum;
}
/**
* The objective is to minimize the distance between all data points
* and all cluster centers. This method returns the summation of all
* these distances
*/
protected function getObjective()
{
$sum = 0.0;
$distance = new Euclidean();
for ($i = 0; $i < $this->clustersNumber; $i++) {
$clust = $this->clusters[$i]->getCoordinates();
for ($k = 0; $k < $this->sampleCount; $k++) {
$point = $this->samples[$k];
$sum += $distance->distance($clust, $point);
}
}
return $sum;
}
/**
* @return array
*/
public function getMembershipMatrix()
{
return $this->membership;
}
/**
* @param array|Point[] $samples
* @return array
*/
public function cluster(array $samples)
{
// Initialize variables, clusters and membership matrix
$this->sampleCount = count($samples);
$this->samples =& $samples;
$this->space = new Space(count($samples[0]));
$this->initClusters();
// Our goal is minimizing the objective value while
// executing the clustering steps at a maximum number of iterations
$lastObjective = 0.0;
$difference = 0.0;
$iterations = 0;
do {
// Update the membership matrix and cluster centers, respectively
$this->updateMembershipMatrix();
$this->updateClusters();
// Calculate the new value of the objective function
$objectiveVal = $this->getObjective();
$difference = abs($lastObjective - $objectiveVal);
$lastObjective = $objectiveVal;
} while ($difference > $this->epsilon && $iterations++ <= $this->maxIterations);
// Attach (hard cluster) each data point to the nearest cluster
for ($k=0; $k<$this->sampleCount; $k++) {
$column = array_column($this->membership, $k);
arsort($column);
reset($column);
$i = key($column);
$cluster = $this->clusters[$i];
$cluster->attach(new Point($this->samples[$k]));
}
// Return grouped samples
$grouped = [];
foreach ($this->clusters as $cluster) {
$grouped[] = $cluster->getPoints();
}
return $grouped;
}
}
@@ -0,0 +1,60 @@
<?php
declare(strict_types=1);
namespace Phpml\Clustering;
use Phpml\Clustering\KMeans\Space;
use Phpml\Exception\InvalidArgumentException;
class KMeans implements Clusterer
{
const INIT_RANDOM = 1;
const INIT_KMEANS_PLUS_PLUS = 2;
/**
* @var int
*/
private $clustersNumber;
/**
* @var int
*/
private $initialization;
/**
* @param int $clustersNumber
* @param int $initialization
*
* @throws InvalidArgumentException
*/
public function __construct(int $clustersNumber, int $initialization = self::INIT_KMEANS_PLUS_PLUS)
{
if ($clustersNumber <= 0) {
throw InvalidArgumentException::invalidClustersNumber();
}
$this->clustersNumber = $clustersNumber;
$this->initialization = $initialization;
}
/**
* @param array $samples
*
* @return array
*/
public function cluster(array $samples)
{
$space = new Space(count($samples[0]));
foreach ($samples as $sample) {
$space->addPoint($sample);
}
$clusters = [];
foreach ($space->cluster($this->clustersNumber, $this->initialization) as $cluster) {
$clusters[] = $cluster->getPoints();
}
return $clusters;
}
}
@@ -0,0 +1,147 @@
<?php
declare(strict_types=1);
namespace Phpml\Clustering\KMeans;
use IteratorAggregate;
use Countable;
use SplObjectStorage;
use LogicException;
class Cluster extends Point implements IteratorAggregate, Countable
{
/**
* @var Space
*/
protected $space;
/**
* @var SplObjectStorage|Point[]
*/
protected $points;
/**
* @param Space $space
* @param array $coordinates
*/
public function __construct(Space $space, array $coordinates)
{
parent::__construct($coordinates);
$this->space = $space;
$this->points = new SplObjectStorage();
}
/**
* @return array
*/
public function getPoints()
{
$points = [];
foreach ($this->points as $point) {
$points[] = $point->toArray();
}
return $points;
}
/**
* @return array
*/
public function toArray()
{
return [
'centroid' => parent::toArray(),
'points' => $this->getPoints(),
];
}
/**
* @param Point $point
*
* @return Point
*
* @throws \LogicException
*/
public function attach(Point $point)
{
if ($point instanceof self) {
throw new LogicException('cannot attach a cluster to another');
}
$this->points->attach($point);
return $point;
}
/**
* @param Point $point
*
* @return Point
*/
public function detach(Point $point)
{
$this->points->detach($point);
return $point;
}
/**
* @param SplObjectStorage $points
*/
public function attachAll(SplObjectStorage $points)
{
$this->points->addAll($points);
}
/**
* @param SplObjectStorage $points
*/
public function detachAll(SplObjectStorage $points)
{
$this->points->removeAll($points);
}
public function updateCentroid()
{
if (!$count = count($this->points)) {
return;
}
$centroid = $this->space->newPoint(array_fill(0, $this->dimension, 0));
foreach ($this->points as $point) {
for ($n = 0; $n < $this->dimension; ++$n) {
$centroid->coordinates[$n] += $point->coordinates[$n];
}
}
for ($n = 0; $n < $this->dimension; ++$n) {
$this->coordinates[$n] = $centroid->coordinates[$n] / $count;
}
}
/**
* @return Point[]|SplObjectStorage
*/
public function getIterator()
{
return $this->points;
}
/**
* @return mixed
*/
public function count()
{
return count($this->points);
}
/**
* @param array $newCoordinates
*/
public function setCoordinates(array $newCoordinates)
{
$this->coordinates = $newCoordinates;
}
}
@@ -0,0 +1,124 @@
<?php
declare(strict_types=1);
namespace Phpml\Clustering\KMeans;
use ArrayAccess;
class Point implements ArrayAccess
{
/**
* @var int
*/
protected $dimension;
/**
* @var array
*/
protected $coordinates;
/**
* @param array $coordinates
*/
public function __construct(array $coordinates)
{
$this->dimension = count($coordinates);
$this->coordinates = $coordinates;
}
/**
* @return array
*/
public function toArray()
{
return $this->coordinates;
}
/**
* @param Point $point
* @param bool $precise
*
* @return int|mixed
*/
public function getDistanceWith(self $point, $precise = true)
{
$distance = 0;
for ($n = 0; $n < $this->dimension; ++$n) {
$difference = $this->coordinates[$n] - $point->coordinates[$n];
$distance += $difference * $difference;
}
return $precise ? sqrt((float) $distance) : $distance;
}
/**
* @param array $points
*
* @return mixed
*/
public function getClosest(array $points)
{
foreach ($points as $point) {
$distance = $this->getDistanceWith($point, false);
if (!isset($minDistance)) {
$minDistance = $distance;
$minPoint = $point;
continue;
}
if ($distance < $minDistance) {
$minDistance = $distance;
$minPoint = $point;
}
}
return $minPoint;
}
/**
* @return array
*/
public function getCoordinates()
{
return $this->coordinates;
}
/**
* @param mixed $offset
*
* @return bool
*/
public function offsetExists($offset)
{
return isset($this->coordinates[$offset]);
}
/**
* @param mixed $offset
*
* @return mixed
*/
public function offsetGet($offset)
{
return $this->coordinates[$offset];
}
/**
* @param mixed $offset
* @param mixed $value
*/
public function offsetSet($offset, $value)
{
$this->coordinates[$offset] = $value;
}
/**
* @param mixed $offset
*/
public function offsetUnset($offset)
{
unset($this->coordinates[$offset]);
}
}
@@ -0,0 +1,259 @@
<?php
declare(strict_types=1);
namespace Phpml\Clustering\KMeans;
use Phpml\Clustering\KMeans;
use SplObjectStorage;
use LogicException;
use InvalidArgumentException;
class Space extends SplObjectStorage
{
/**
* @var int
*/
protected $dimension;
/**
* @param $dimension
*/
public function __construct($dimension)
{
if ($dimension < 1) {
throw new LogicException('a space dimension cannot be null or negative');
}
$this->dimension = $dimension;
}
/**
* @return array
*/
public function toArray()
{
$points = [];
foreach ($this as $point) {
$points[] = $point->toArray();
}
return ['points' => $points];
}
/**
* @param array $coordinates
*
* @return Point
*/
public function newPoint(array $coordinates)
{
if (count($coordinates) != $this->dimension) {
throw new LogicException('('.implode(',', $coordinates).') is not a point of this space');
}
return new Point($coordinates);
}
/**
* @param array $coordinates
* @param null $data
*/
public function addPoint(array $coordinates, $data = null)
{
$this->attach($this->newPoint($coordinates), $data);
}
/**
* @param Point $point
* @param null $data
*/
public function attach($point, $data = null)
{
if (!$point instanceof Point) {
throw new InvalidArgumentException('can only attach points to spaces');
}
parent::attach($point, $data);
}
/**
* @return int
*/
public function getDimension()
{
return $this->dimension;
}
/**
* @return array|bool
*/
public function getBoundaries()
{
if (!count($this)) {
return false;
}
$min = $this->newPoint(array_fill(0, $this->dimension, null));
$max = $this->newPoint(array_fill(0, $this->dimension, null));
foreach ($this as $point) {
for ($n = 0; $n < $this->dimension; ++$n) {
($min[$n] > $point[$n] || $min[$n] === null) && $min[$n] = $point[$n];
($max[$n] < $point[$n] || $max[$n] === null) && $max[$n] = $point[$n];
}
}
return [$min, $max];
}
/**
* @param Point $min
* @param Point $max
*
* @return Point
*/
public function getRandomPoint(Point $min, Point $max)
{
$point = $this->newPoint(array_fill(0, $this->dimension, null));
for ($n = 0; $n < $this->dimension; ++$n) {
$point[$n] = random_int($min[$n], $max[$n]);
}
return $point;
}
/**
* @param int $clustersNumber
* @param int $initMethod
*
* @return array|Cluster[]
*/
public function cluster(int $clustersNumber, int $initMethod = KMeans::INIT_RANDOM)
{
$clusters = $this->initializeClusters($clustersNumber, $initMethod);
do {
} while (!$this->iterate($clusters));
return $clusters;
}
/**
* @param $clustersNumber
* @param $initMethod
*
* @return array|Cluster[]
*/
protected function initializeClusters(int $clustersNumber, int $initMethod)
{
switch ($initMethod) {
case KMeans::INIT_RANDOM:
$clusters = $this->initializeRandomClusters($clustersNumber);
break;
case KMeans::INIT_KMEANS_PLUS_PLUS:
$clusters = $this->initializeKMPPClusters($clustersNumber);
break;
}
$clusters[0]->attachAll($this);
return $clusters;
}
/**
* @param $clusters
*
* @return bool
*/
protected function iterate($clusters)
{
$convergence = true;
$attach = new SplObjectStorage();
$detach = new SplObjectStorage();
foreach ($clusters as $cluster) {
foreach ($cluster as $point) {
$closest = $point->getClosest($clusters);
if ($closest !== $cluster) {
isset($attach[$closest]) || $attach[$closest] = new SplObjectStorage();
isset($detach[$cluster]) || $detach[$cluster] = new SplObjectStorage();
$attach[$closest]->attach($point);
$detach[$cluster]->attach($point);
$convergence = false;
}
}
}
foreach ($attach as $cluster) {
$cluster->attachAll($attach[$cluster]);
}
foreach ($detach as $cluster) {
$cluster->detachAll($detach[$cluster]);
}
foreach ($clusters as $cluster) {
$cluster->updateCentroid();
}
return $convergence;
}
/**
* @param int $clustersNumber
*
* @return array
*/
private function initializeRandomClusters(int $clustersNumber)
{
$clusters = [];
list($min, $max) = $this->getBoundaries();
for ($n = 0; $n < $clustersNumber; ++$n) {
$clusters[] = new Cluster($this, $this->getRandomPoint($min, $max)->getCoordinates());
}
return $clusters;
}
/**
* @param int $clustersNumber
*
* @return array
*/
protected function initializeKMPPClusters(int $clustersNumber)
{
$clusters = [];
$this->rewind();
$clusters[] = new Cluster($this, $this->current()->getCoordinates());
$distances = new SplObjectStorage();
for ($i = 1; $i < $clustersNumber; ++$i) {
$sum = 0;
foreach ($this as $point) {
$distance = $point->getDistanceWith($point->getClosest($clusters));
$sum += $distances[$point] = $distance;
}
$sum = random_int(0, (int) $sum);
foreach ($this as $point) {
if (($sum -= $distances[$point]) > 0) {
continue;
}
$clusters[] = new Cluster($this, $point->getCoordinates());
break;
}
}
return $clusters;
}
}
@@ -0,0 +1,30 @@
<?php
declare(strict_types=1);
namespace Phpml\CrossValidation;
use Phpml\Dataset\Dataset;
class RandomSplit extends Split
{
/**
* @param Dataset $dataset
* @param float $testSize
*/
protected function splitDataset(Dataset $dataset, float $testSize)
{
$samples = $dataset->getSamples();
$labels = $dataset->getTargets();
$datasetSize = count($samples);
$testCount = count($this->testSamples);
for ($i = $datasetSize; $i > 0; --$i) {
$key = mt_rand(0, $datasetSize - 1);
$setName = (count($this->testSamples) - $testCount) / $datasetSize >= $testSize ? 'train' : 'test';
$this->{$setName.'Samples'}[] = $samples[$key];
$this->{$setName.'Labels'}[] = $labels[$key];
}
}
}
@@ -0,0 +1,94 @@
<?php
declare(strict_types=1);
namespace Phpml\CrossValidation;
use Phpml\Dataset\Dataset;
use Phpml\Exception\InvalidArgumentException;
abstract class Split
{
/**
* @var array
*/
protected $trainSamples = [];
/**
* @var array
*/
protected $testSamples = [];
/**
* @var array
*/
protected $trainLabels = [];
/**
* @var array
*/
protected $testLabels = [];
/**
* @param Dataset $dataset
* @param float $testSize
* @param int $seed
*
* @throws InvalidArgumentException
*/
public function __construct(Dataset $dataset, float $testSize = 0.3, int $seed = null)
{
if (0 >= $testSize || 1 <= $testSize) {
throw InvalidArgumentException::percentNotInRange('testSize');
}
$this->seedGenerator($seed);
$this->splitDataset($dataset, $testSize);
}
abstract protected function splitDataset(Dataset $dataset, float $testSize);
/**
* @return array
*/
public function getTrainSamples()
{
return $this->trainSamples;
}
/**
* @return array
*/
public function getTestSamples()
{
return $this->testSamples;
}
/**
* @return array
*/
public function getTrainLabels()
{
return $this->trainLabels;
}
/**
* @return array
*/
public function getTestLabels()
{
return $this->testLabels;
}
/**
* @param int|null $seed
*/
protected function seedGenerator(int $seed = null)
{
if (null === $seed) {
mt_srand();
} else {
mt_srand($seed);
}
}
}
@@ -0,0 +1,62 @@
<?php
declare(strict_types=1);
namespace Phpml\CrossValidation;
use Phpml\Dataset\ArrayDataset;
use Phpml\Dataset\Dataset;
class StratifiedRandomSplit extends RandomSplit
{
/**
* @param Dataset $dataset
* @param float $testSize
*/
protected function splitDataset(Dataset $dataset, float $testSize)
{
$datasets = $this->splitByTarget($dataset);
foreach ($datasets as $targetSet) {
parent::splitDataset($targetSet, $testSize);
}
}
/**
* @param Dataset $dataset
*
* @return Dataset[]|array
*/
private function splitByTarget(Dataset $dataset): array
{
$targets = $dataset->getTargets();
$samples = $dataset->getSamples();
$uniqueTargets = array_unique($targets);
$split = array_combine($uniqueTargets, array_fill(0, count($uniqueTargets), []));
foreach ($samples as $key => $sample) {
$split[$targets[$key]][] = $sample;
}
$datasets = $this->createDatasets($uniqueTargets, $split);
return $datasets;
}
/**
* @param array $uniqueTargets
* @param array $split
*
* @return array
*/
private function createDatasets(array $uniqueTargets, array $split): array
{
$datasets = [];
foreach ($uniqueTargets as $target) {
$datasets[$target] = new ArrayDataset($split[$target], array_fill(0, count($split[$target]), $target));
}
return $datasets;
}
}
@@ -0,0 +1,52 @@
<?php
declare(strict_types=1);
namespace Phpml\Dataset;
use Phpml\Exception\InvalidArgumentException;
class ArrayDataset implements Dataset
{
/**
* @var array
*/
protected $samples = [];
/**
* @var array
*/
protected $targets = [];
/**
* @param array $samples
* @param array $targets
*
* @throws InvalidArgumentException
*/
public function __construct(array $samples, array $targets)
{
if (count($samples) != count($targets)) {
throw InvalidArgumentException::arraySizeNotMatch();
}
$this->samples = $samples;
$this->targets = $targets;
}
/**
* @return array
*/
public function getSamples(): array
{
return $this->samples;
}
/**
* @return array
*/
public function getTargets(): array
{
return $this->targets;
}
}
@@ -0,0 +1,54 @@
<?php
declare(strict_types=1);
namespace Phpml\Dataset;
use Phpml\Exception\FileException;
class CsvDataset extends ArrayDataset
{
/**
* @var array
*/
protected $columnNames;
/**
* @param string $filepath
* @param int $features
* @param bool $headingRow
*
* @throws FileException
*/
public function __construct(string $filepath, int $features, bool $headingRow = true, string $delimiter = ',')
{
if (!file_exists($filepath)) {
throw FileException::missingFile(basename($filepath));
}
if (false === $handle = fopen($filepath, 'rb')) {
throw FileException::cantOpenFile(basename($filepath));
}
if ($headingRow) {
$data = fgetcsv($handle, 1000, $delimiter);
$this->columnNames = array_slice($data, 0, $features);
} else {
$this->columnNames = range(0, $features - 1);
}
while (($data = fgetcsv($handle, 1000, $delimiter)) !== false) {
$this->samples[] = array_slice($data, 0, $features);
$this->targets[] = $data[$features];
}
fclose($handle);
}
/**
* @return array
*/
public function getColumnNames()
{
return $this->columnNames;
}
}
@@ -0,0 +1,18 @@
<?php
declare(strict_types=1);
namespace Phpml\Dataset;
interface Dataset
{
/**
* @return array
*/
public function getSamples(): array;
/**
* @return array
*/
public function getTargets(): array;
}
@@ -0,0 +1,28 @@
<?php
declare(strict_types=1);
namespace Phpml\Dataset\Demo;
use Phpml\Dataset\CsvDataset;
/**
* Classes: 6
* Samples per class:
* 70 float processed building windows
* 17 float processed vehicle windows
* 76 non-float processed building windows
* 13 containers
* 9 tableware
* 29 headlamps
* Samples total: 214
* Features per sample: 9.
*/
class GlassDataset extends CsvDataset
{
public function __construct()
{
$filepath = __DIR__.'/../../../../data/glass.csv';
parent::__construct($filepath, 9, true);
}
}
@@ -0,0 +1,22 @@
<?php
declare(strict_types=1);
namespace Phpml\Dataset\Demo;
use Phpml\Dataset\CsvDataset;
/**
* Classes: 3
* Samples per class: 50
* Samples total: 150
* Features per sample: 4.
*/
class IrisDataset extends CsvDataset
{
public function __construct()
{
$filepath = __DIR__.'/../../../../data/iris.csv';
parent::__construct($filepath, 4, true);
}
}
@@ -0,0 +1,22 @@
<?php
declare(strict_types=1);
namespace Phpml\Dataset\Demo;
use Phpml\Dataset\CsvDataset;
/**
* Classes: 3
* Samples per class: class 1 59; class 2 71; class 3 48
* Samples total: 178
* Features per sample: 13.
*/
class WineDataset extends CsvDataset
{
public function __construct()
{
$filepath = __DIR__.'/../../../../data/wine.csv';
parent::__construct($filepath, 13, true);
}
}
@@ -0,0 +1,47 @@
<?php
declare(strict_types=1);
namespace Phpml\Dataset;
use Phpml\Exception\DatasetException;
class FilesDataset extends ArrayDataset
{
/**
* @param string $rootPath
*
* @throws DatasetException
*/
public function __construct(string $rootPath)
{
if (!is_dir($rootPath)) {
throw DatasetException::missingFolder($rootPath);
}
$this->scanRootPath($rootPath);
}
/**
* @param string $rootPath
*/
private function scanRootPath(string $rootPath)
{
foreach (glob($rootPath.DIRECTORY_SEPARATOR.'*', GLOB_ONLYDIR) as $dir) {
$this->scanDir($dir);
}
}
/**
* @param string $dir
*/
private function scanDir(string $dir)
{
$target = basename($dir);
foreach (array_filter(glob($dir.DIRECTORY_SEPARATOR.'*'), 'is_file') as $file) {
$this->samples[] = [file_get_contents($file)];
$this->targets[] = $target;
}
}
}
@@ -0,0 +1,246 @@
<?php
declare(strict_types=1);
namespace Phpml\DimensionReduction;
use Phpml\Math\Distance\Euclidean;
use Phpml\Math\Distance\Manhattan;
use Phpml\Math\Matrix;
class KernelPCA extends PCA
{
const KERNEL_RBF = 1;
const KERNEL_SIGMOID = 2;
const KERNEL_LAPLACIAN = 3;
const KERNEL_LINEAR = 4;
/**
* Selected kernel function
*
* @var int
*/
protected $kernel;
/**
* Gamma value used by the kernel
*
* @var float
*/
protected $gamma;
/**
* Original dataset used to fit KernelPCA
*
* @var array
*/
protected $data;
/**
* Kernel principal component analysis (KernelPCA) is an extension of PCA using
* techniques of kernel methods. It is more suitable for data that involves
* vectors that are not linearly separable<br><br>
* Example: <b>$kpca = new KernelPCA(KernelPCA::KERNEL_RBF, null, 2, 15.0);</b>
* will initialize the algorithm with an RBF kernel having the gamma parameter as 15,0. <br>
* This transformation will return the same number of rows with only <i>2</i> columns.
*
* @param int $kernel
* @param float $totalVariance Total variance to be preserved if numFeatures is not given
* @param int $numFeatures Number of columns to be returned
* @param float $gamma Gamma parameter is used with RBF and Sigmoid kernels
*
* @throws \Exception
*/
public function __construct(int $kernel = self::KERNEL_RBF, $totalVariance = null, $numFeatures = null, $gamma = null)
{
$availableKernels = [self::KERNEL_RBF, self::KERNEL_SIGMOID, self::KERNEL_LAPLACIAN, self::KERNEL_LINEAR];
if (! in_array($kernel, $availableKernels)) {
throw new \Exception("KernelPCA can be initialized with the following kernels only: Linear, RBF, Sigmoid and Laplacian");
}
parent::__construct($totalVariance, $numFeatures);
$this->kernel = $kernel;
$this->gamma = $gamma;
}
/**
* Takes a data and returns a lower dimensional version
* of this data while preserving $totalVariance or $numFeatures. <br>
* $data is an n-by-m matrix and returned array is
* n-by-k matrix where k <= m
*
* @param array $data
*
* @return array
*/
public function fit(array $data)
{
$numRows = count($data);
$this->data = $data;
if ($this->gamma === null) {
$this->gamma = 1.0 / $numRows;
}
$matrix = $this->calculateKernelMatrix($this->data, $numRows);
$matrix = $this->centerMatrix($matrix, $numRows);
list($this->eigValues, $this->eigVectors) = $this->eigenDecomposition($matrix, $numRows);
$this->fit = true;
return Matrix::transposeArray($this->eigVectors);
}
/**
* Calculates similarity matrix by use of selected kernel function<br>
* An n-by-m matrix is given and an n-by-n matrix is returned
*
* @param array $data
* @param int $numRows
*
* @return array
*/
protected function calculateKernelMatrix(array $data, int $numRows)
{
$kernelFunc = $this->getKernel();
$matrix = [];
for ($i=0; $i < $numRows; $i++) {
for ($k=0; $k < $numRows; $k++) {
if ($i <= $k) {
$matrix[$i][$k] = $kernelFunc($data[$i], $data[$k]);
} else {
$matrix[$i][$k] = $matrix[$k][$i];
}
}
}
return $matrix;
}
/**
* Kernel matrix is centered in its original space by using the following
* conversion:
*
* K = K N.K K.N + N.K.N where N is n-by-n matrix filled with 1/n
*
* @param array $matrix
* @param int $n
*/
protected function centerMatrix(array $matrix, int $n)
{
$N = array_fill(0, $n, array_fill(0, $n, 1.0/$n));
$N = new Matrix($N, false);
$K = new Matrix($matrix, false);
// K.N (This term is repeated so we cache it once)
$K_N = $K->multiply($N);
// N.K
$N_K = $N->multiply($K);
// N.K.N
$N_K_N = $N->multiply($K_N);
return $K->subtract($N_K)
->subtract($K_N)
->add($N_K_N)
->toArray();
}
/**
* Returns the callable kernel function
*
* @return \Closure
*/
protected function getKernel()
{
switch ($this->kernel) {
case self::KERNEL_LINEAR:
// k(x,y) = xT.y
return function ($x, $y) {
return Matrix::dot($x, $y)[0];
};
case self::KERNEL_RBF:
// k(x,y)=exp(-γ.|x-y|) where |..| is Euclidean distance
$dist = new Euclidean();
return function ($x, $y) use ($dist) {
return exp(-$this->gamma * $dist->sqDistance($x, $y));
};
case self::KERNEL_SIGMOID:
// k(x,y)=tanh(γ.xT.y+c0) where c0=1
return function ($x, $y) {
$res = Matrix::dot($x, $y)[0] + 1.0;
return tanh($this->gamma * $res);
};
case self::KERNEL_LAPLACIAN:
// k(x,y)=exp(-γ.|x-y|) where |..| is Manhattan distance
$dist = new Manhattan();
return function ($x, $y) use ($dist) {
return exp(-$this->gamma * $dist->distance($x, $y));
};
}
}
/**
* @param array $sample
*
* @return array
*/
protected function getDistancePairs(array $sample)
{
$kernel = $this->getKernel();
$pairs = [];
foreach ($this->data as $row) {
$pairs[] = $kernel($row, $sample);
}
return $pairs;
}
/**
* @param array $pairs
*
* @return array
*/
protected function projectSample(array $pairs)
{
// Normalize eigenvectors by eig = eigVectors / eigValues
$func = function ($eigVal, $eigVect) {
$m = new Matrix($eigVect, false);
$a = $m->divideByScalar($eigVal)->toArray();
return $a[0];
};
$eig = array_map($func, $this->eigValues, $this->eigVectors);
// return k.dot(eig)
return Matrix::dot($pairs, $eig);
}
/**
* Transforms the given sample to a lower dimensional vector by using
* the variables obtained during the last run of <code>fit</code>.
*
* @param array $sample
*
* @return array
*/
public function transform(array $sample)
{
if (!$this->fit) {
throw new \Exception("KernelPCA has not been fitted with respect to original dataset, please run KernelPCA::fit() first");
}
if (is_array($sample[0])) {
throw new \Exception("KernelPCA::transform() accepts only one-dimensional arrays");
}
$pairs = $this->getDistancePairs($sample);
return $this->projectSample($pairs);
}
}
@@ -0,0 +1,228 @@
<?php
declare(strict_types=1);
namespace Phpml\DimensionReduction;
use Phpml\Math\LinearAlgebra\EigenvalueDecomposition;
use Phpml\Math\Statistic\Covariance;
use Phpml\Math\Statistic\Mean;
use Phpml\Math\Matrix;
class PCA
{
/**
* Total variance to be conserved after the reduction
*
* @var float
*/
public $totalVariance = 0.9;
/**
* Number of features to be preserved after the reduction
*
* @var int
*/
public $numFeatures = null;
/**
* Temporary storage for mean values for each dimension in given data
*
* @var array
*/
protected $means = [];
/**
* Eigenvectors of the covariance matrix
*
* @var array
*/
protected $eigVectors = [];
/**
* Top eigenValues of the covariance matrix
*
* @var type
*/
protected $eigValues = [];
/**
* @var bool
*/
protected $fit = false;
/**
* PCA (Principal Component Analysis) used to explain given
* data with lower number of dimensions. This analysis transforms the
* data to a lower dimensional version of it by conserving a proportion of total variance
* within the data. It is a lossy data compression technique.<br>
*
* @param float $totalVariance Total explained variance to be preserved
* @param int $numFeatures Number of features to be preserved
*
* @throws \Exception
*/
public function __construct($totalVariance = null, $numFeatures = null)
{
if ($totalVariance !== null && ($totalVariance < 0.1 || $totalVariance > 0.99)) {
throw new \Exception("Total variance can be a value between 0.1 and 0.99");
}
if ($numFeatures !== null && $numFeatures <= 0) {
throw new \Exception("Number of features to be preserved should be greater than 0");
}
if ($totalVariance !== null && $numFeatures !== null) {
throw new \Exception("Either totalVariance or numFeatures should be specified in order to run the algorithm");
}
if ($numFeatures !== null) {
$this->numFeatures = $numFeatures;
}
if ($totalVariance !== null) {
$this->totalVariance = $totalVariance;
}
}
/**
* Takes a data and returns a lower dimensional version
* of this data while preserving $totalVariance or $numFeatures. <br>
* $data is an n-by-m matrix and returned array is
* n-by-k matrix where k <= m
*
* @param array $data
*
* @return array
*/
public function fit(array $data)
{
$n = count($data[0]);
$data = $this->normalize($data, $n);
$covMatrix = Covariance::covarianceMatrix($data, array_fill(0, $n, 0));
list($this->eigValues, $this->eigVectors) = $this->eigenDecomposition($covMatrix, $n);
$this->fit = true;
return $this->reduce($data);
}
/**
* @param array $data
* @param int $n
*/
protected function calculateMeans(array $data, int $n)
{
// Calculate means for each dimension
$this->means = [];
for ($i=0; $i < $n; $i++) {
$column = array_column($data, $i);
$this->means[] = Mean::arithmetic($column);
}
}
/**
* Normalization of the data includes subtracting mean from
* each dimension therefore dimensions will be centered to zero
*
* @param array $data
* @param int $n
*
* @return array
*/
protected function normalize(array $data, int $n)
{
if (empty($this->means)) {
$this->calculateMeans($data, $n);
}
// Normalize data
foreach ($data as $i => $row) {
for ($k=0; $k < $n; $k++) {
$data[$i][$k] -= $this->means[$k];
}
}
return $data;
}
/**
* Calculates eigenValues and eigenVectors of the given matrix. Returns
* top eigenVectors along with the largest eigenValues. The total explained variance
* of these eigenVectors will be no less than desired $totalVariance value
*
* @param array $matrix
* @param int $n
*
* @return array
*/
protected function eigenDecomposition(array $matrix, int $n)
{
$eig = new EigenvalueDecomposition($matrix);
$eigVals = $eig->getRealEigenvalues();
$eigVects= $eig->getEigenvectors();
$totalEigVal = array_sum($eigVals);
// Sort eigenvalues in descending order
arsort($eigVals);
$explainedVar = 0.0;
$vectors = [];
$values = [];
foreach ($eigVals as $i => $eigVal) {
$explainedVar += $eigVal / $totalEigVal;
$vectors[] = $eigVects[$i];
$values[] = $eigVal;
if ($this->numFeatures !== null) {
if (count($vectors) == $this->numFeatures) {
break;
}
} else {
if ($explainedVar >= $this->totalVariance) {
break;
}
}
}
return [$values, $vectors];
}
/**
* Returns the reduced data
*
* @param array $data
*
* @return array
*/
protected function reduce(array $data)
{
$m1 = new Matrix($data);
$m2 = new Matrix($this->eigVectors);
return $m1->multiply($m2->transpose())->toArray();
}
/**
* Transforms the given sample to a lower dimensional vector by using
* the eigenVectors obtained in the last run of <code>fit</code>.
*
* @param array $sample
*
* @return array
*/
public function transform(array $sample)
{
if (!$this->fit) {
throw new \Exception("PCA has not been fitted with respect to original dataset, please run PCA::fit() first");
}
if (! is_array($sample[0])) {
$sample = [$sample];
}
$sample = $this->normalize($sample, count($sample[0]));
return $this->reduce($sample);
}
}
@@ -0,0 +1,21 @@
<?php
declare(strict_types=1);
namespace Phpml;
interface Estimator
{
/**
* @param array $samples
* @param array $targets
*/
public function train(array $samples, array $targets);
/**
* @param array $samples
*
* @return mixed
*/
public function predict(array $samples);
}
@@ -0,0 +1,19 @@
<?php
declare(strict_types=1);
namespace Phpml\Exception;
class DatasetException extends \Exception
{
/**
* @param string $path
*
* @return DatasetException
*/
public static function missingFolder(string $path)
{
return new self(sprintf('Dataset root folder "%s" missing.', $path));
}
}
@@ -0,0 +1,39 @@
<?php
declare(strict_types=1);
namespace Phpml\Exception;
class FileException extends \Exception
{
/**
* @param string $filepath
*
* @return FileException
*/
public static function missingFile(string $filepath)
{
return new self(sprintf('File "%s" missing.', $filepath));
}
/**
* @param string $filepath
*
* @return FileException
*/
public static function cantOpenFile(string $filepath)
{
return new self(sprintf('File "%s" can\'t be open.', $filepath));
}
/**
* @param string $filepath
*
* @return FileException
*/
public static function cantSaveFile(string $filepath)
{
return new self(sprintf('File "%s" can\'t be saved.', $filepath));
}
}
@@ -0,0 +1,94 @@
<?php
declare(strict_types=1);
namespace Phpml\Exception;
class InvalidArgumentException extends \Exception
{
/**
* @return InvalidArgumentException
*/
public static function arraySizeNotMatch()
{
return new self('Size of given arrays not match');
}
/**
* @param $name
*
* @return InvalidArgumentException
*/
public static function percentNotInRange($name)
{
return new self(sprintf('%s must be between 0.0 and 1.0', $name));
}
/**
* @return InvalidArgumentException
*/
public static function arrayCantBeEmpty()
{
return new self('The array has zero elements');
}
/**
* @param int $minimumSize
*
* @return InvalidArgumentException
*/
public static function arraySizeToSmall($minimumSize = 2)
{
return new self(sprintf('The array must have at least %s elements', $minimumSize));
}
/**
* @return InvalidArgumentException
*/
public static function matrixDimensionsDidNotMatch()
{
return new self('Matrix dimensions did not match');
}
/**
* @return InvalidArgumentException
*/
public static function inconsistentMatrixSupplied()
{
return new self('Inconsistent matrix applied');
}
/**
* @return InvalidArgumentException
*/
public static function invalidClustersNumber()
{
return new self('Invalid clusters number');
}
/**
* @param string $language
*
* @return InvalidArgumentException
*/
public static function invalidStopWordsLanguage(string $language)
{
return new self(sprintf('Can\'t find %s language for StopWords', $language));
}
/**
* @return InvalidArgumentException
*/
public static function invalidLayerNodeClass()
{
return new self('Layer node class must implement Node interface');
}
/**
* @return InvalidArgumentException
*/
public static function invalidLayersNumber()
{
return new self('Provide at least 2 layers: 1 input and 1 output');
}
}
@@ -0,0 +1,32 @@
<?php
declare(strict_types=1);
namespace Phpml\Exception;
class MatrixException extends \Exception
{
/**
* @return MatrixException
*/
public static function notSquareMatrix()
{
return new self('Matrix is not square matrix');
}
/**
* @return MatrixException
*/
public static function columnOutOfRange()
{
return new self('Column out of range');
}
/**
* @return MatrixException
*/
public static function singularMatrix()
{
return new self('Matrix is singular');
}
}
@@ -0,0 +1,16 @@
<?php
declare(strict_types=1);
namespace Phpml\Exception;
class NormalizerException extends \Exception
{
/**
* @return NormalizerException
*/
public static function unknownNorm()
{
return new self('Unknown norm supplied.');
}
}
@@ -0,0 +1,29 @@
<?php
declare(strict_types=1);
namespace Phpml\Exception;
class SerializeException extends \Exception
{
/**
* @param string $filepath
*
* @return SerializeException
*/
public static function cantUnserialize(string $filepath)
{
return new self(sprintf('"%s" can not be unserialized.', $filepath));
}
/**
* @param string $classname
*
* @return SerializeException
*/
public static function cantSerialize(string $classname)
{
return new self(sprintf('Class "%s" can not be serialized.', $classname));
}
}
@@ -0,0 +1,51 @@
<?php
declare(strict_types=1);
namespace Phpml\FeatureExtraction;
use Phpml\Exception\InvalidArgumentException;
class StopWords
{
/**
* @var array
*/
protected $stopWords;
/**
* @param array $stopWords
*/
public function __construct(array $stopWords)
{
$this->stopWords = array_fill_keys($stopWords, true);
}
/**
* @param string $token
*
* @return bool
*/
public function isStopWord(string $token): bool
{
return isset($this->stopWords[$token]);
}
/**
* @param string $language
*
* @return StopWords
*
* @throws InvalidArgumentException
*/
public static function factory($language = 'English'): StopWords
{
$className = __NAMESPACE__."\\StopWords\\$language";
if (!class_exists($className)) {
throw InvalidArgumentException::invalidStopWordsLanguage($language);
}
return new $className();
}
}
@@ -0,0 +1,33 @@
<?php
declare(strict_types=1);
namespace Phpml\FeatureExtraction\StopWords;
use Phpml\FeatureExtraction\StopWords;
final class English extends StopWords
{
/**
* @var array
*/
protected $stopWords = [
'a', 'about', 'above', 'after', 'again', 'against', 'all', 'am', 'an', 'and', 'any', 'are', 'aren\'t', 'as', 'at', 'be', 'because',
'been', 'before', 'being', 'below', 'between', 'both', 'but', 'by', 'can\'t', 'cannot', 'could', 'couldn\'t', 'did', 'didn\'t',
'do', 'does', 'doesn\'t', 'doing', 'don\'t', 'down', 'during', 'each', 'few', 'for', 'from', 'further', 'had', 'hadn\'t', 'has',
'hasn\'t', 'have', 'haven\'t', 'having', 'he', 'he\'d', 'he\'ll', 'he\'s', 'her', 'here', 'here\'s', 'hers', 'herself', 'him',
'himself', 'his', 'how', 'how\'s', 'i', 'i\'d', 'i\'ll', 'i\'m', 'i\'ve', 'if', 'in', 'into', 'is', 'isn\'t', 'it', 'it\'s', 'its',
'itself', 'let\'s', 'me', 'more', 'most', 'mustn\'t', 'my', 'myself', 'no', 'nor', 'not', 'of', 'off', 'on', 'once', 'only', 'or',
'other', 'ought', 'our', 'oursourselves', 'out', 'over', 'own', 'same', 'shan\'t', 'she', 'she\'d', 'she\'ll', 'she\'s', 'should',
'shouldn\'t', 'so', 'some', 'such', 'than', 'that', 'that\'s', 'the', 'their', 'theirs', 'them', 'themselves', 'then', 'there',
'there\'s', 'these', 'they', 'they\'d', 'they\'ll', 'they\'re', 'they\'ve', 'this', 'those', 'through', 'to', 'too', 'under',
'until', 'up', 'very', 'was', 'wasn\'t', 'we', 'we\'d', 'we\'ll', 'we\'re', 'we\'ve', 'were', 'weren\'t', 'what', 'what\'s',
'when', 'when\'s', 'where', 'where\'s', 'which', 'while', 'who', 'who\'s', 'whom', 'why', 'why\'s', 'with', 'won\'t', 'would',
'wouldn\'t', 'you', 'you\'d', 'you\'ll', 'you\'re', 'you\'ve', 'your', 'yours', 'yourself', 'yourselves',
];
public function __construct()
{
parent::__construct($this->stopWords);
}
}
@@ -0,0 +1,30 @@
<?php
declare(strict_types=1);
namespace Phpml\FeatureExtraction\StopWords;
use Phpml\FeatureExtraction\StopWords;
final class Polish extends StopWords
{
/**
* @var array
*/
protected $stopWords = [
'ach', 'aj', 'albo', 'bardzo', 'bez', 'bo', 'być', 'ci', 'cię', 'ciebie', 'co', 'czy', 'daleko', 'dla', 'dlaczego', 'dlatego',
'do', 'dobrze', 'dokąd', 'dość', 'dużo', 'dwa', 'dwaj', 'dwie', 'dwoje', 'dziś', 'dzisiaj', 'gdyby', 'gdzie', 'go', 'ich', 'ile',
'im', 'inny', 'ja', 'ją', 'jak', 'jakby', 'jaki', 'je', 'jeden', 'jedna', 'jedno', 'jego', 'jej', 'jemu', 'jeśli', 'jest', 'jestem',
'jeżeli', 'już', 'każdy', 'kiedy', 'kierunku', 'kto', 'ku', 'lub', 'ma', 'mają', 'mam', 'mi', 'mną', 'mnie', 'moi', 'mój', 'moja',
'moje', 'może', 'mu', 'my', 'na', 'nam', 'nami', 'nas', 'nasi', 'nasz', 'nasza', 'nasze', 'natychmiast', 'nią', 'nic', 'nich',
'nie', 'niego', 'niej', 'niemu', 'nigdy', 'nim', 'nimi', 'niż', 'obok', 'od', 'około', 'on', 'ona', 'one', 'oni', 'ono', 'owszem',
'po', 'pod', 'ponieważ', 'przed', 'przedtem', 'są', 'sam', 'sama', 'się', 'skąd', 'tak', 'taki', 'tam', 'ten', 'to', 'tobą', 'tobie',
'tu', 'tutaj', 'twoi', 'twój', 'twoja', 'twoje', 'ty', 'wam', 'wami', 'was', 'wasi', 'wasz', 'wasza', 'wasze', 'we', 'więc',
'wszystko', 'wtedy', 'wy', 'żaden', 'zawsze', 'że',
];
public function __construct()
{
parent::__construct($this->stopWords);
}
}
@@ -0,0 +1,66 @@
<?php
declare(strict_types=1);
namespace Phpml\FeatureExtraction;
use Phpml\Transformer;
class TfIdfTransformer implements Transformer
{
/**
* @var array
*/
private $idf;
/**
* @param array $samples
*/
public function __construct(array $samples = null)
{
if ($samples) {
$this->fit($samples);
}
}
/**
* @param array $samples
*/
public function fit(array $samples)
{
$this->countTokensFrequency($samples);
$count = count($samples);
foreach ($this->idf as &$value) {
$value = log((float)($count / $value), 10.0);
}
}
/**
* @param array $samples
*/
public function transform(array &$samples)
{
foreach ($samples as &$sample) {
foreach ($sample as $index => &$feature) {
$feature *= $this->idf[$index];
}
}
}
/**
* @param array $samples
*/
private function countTokensFrequency(array $samples)
{
$this->idf = array_fill_keys(array_keys($samples[0]), 0);
foreach ($samples as $sample) {
foreach ($sample as $index => $count) {
if ($count > 0) {
++$this->idf[$index];
}
}
}
}
}
@@ -0,0 +1,214 @@
<?php
declare(strict_types=1);
namespace Phpml\FeatureExtraction;
use Phpml\Tokenization\Tokenizer;
use Phpml\Transformer;
class TokenCountVectorizer implements Transformer
{
/**
* @var Tokenizer
*/
private $tokenizer;
/**
* @var StopWords
*/
private $stopWords;
/**
* @var float
*/
private $minDF;
/**
* @var array
*/
private $vocabulary;
/**
* @var array
*/
private $frequencies;
/**
* @param Tokenizer $tokenizer
* @param StopWords $stopWords
* @param float $minDF
*/
public function __construct(Tokenizer $tokenizer, StopWords $stopWords = null, float $minDF = 0.0)
{
$this->tokenizer = $tokenizer;
$this->stopWords = $stopWords;
$this->minDF = $minDF;
$this->vocabulary = [];
$this->frequencies = [];
}
/**
* @param array $samples
*/
public function fit(array $samples)
{
$this->buildVocabulary($samples);
}
/**
* @param array $samples
*/
public function transform(array &$samples)
{
foreach ($samples as &$sample) {
$this->transformSample($sample);
}
$this->checkDocumentFrequency($samples);
}
/**
* @return array
*/
public function getVocabulary()
{
return array_flip($this->vocabulary);
}
/**
* @param array $samples
*/
private function buildVocabulary(array &$samples)
{
foreach ($samples as $index => $sample) {
$tokens = $this->tokenizer->tokenize($sample);
foreach ($tokens as $token) {
$this->addTokenToVocabulary($token);
}
}
}
/**
* @param string $sample
*/
private function transformSample(string &$sample)
{
$counts = [];
$tokens = $this->tokenizer->tokenize($sample);
foreach ($tokens as $token) {
$index = $this->getTokenIndex($token);
if (false !== $index) {
$this->updateFrequency($token);
if (!isset($counts[$index])) {
$counts[$index] = 0;
}
++$counts[$index];
}
}
foreach ($this->vocabulary as $index) {
if (!isset($counts[$index])) {
$counts[$index] = 0;
}
}
ksort($counts);
$sample = $counts;
}
/**
* @param string $token
*
* @return int|bool
*/
private function getTokenIndex(string $token)
{
if ($this->isStopWord($token)) {
return false;
}
return $this->vocabulary[$token] ?? false;
}
/**
* @param string $token
*/
private function addTokenToVocabulary(string $token)
{
if ($this->isStopWord($token)) {
return;
}
if (!isset($this->vocabulary[$token])) {
$this->vocabulary[$token] = count($this->vocabulary);
}
}
/**
* @param string $token
*
* @return bool
*/
private function isStopWord(string $token): bool
{
return $this->stopWords && $this->stopWords->isStopWord($token);
}
/**
* @param string $token
*/
private function updateFrequency(string $token)
{
if (!isset($this->frequencies[$token])) {
$this->frequencies[$token] = 0;
}
++$this->frequencies[$token];
}
/**
* @param array $samples
*/
private function checkDocumentFrequency(array &$samples)
{
if ($this->minDF > 0) {
$beyondMinimum = $this->getBeyondMinimumIndexes(count($samples));
foreach ($samples as &$sample) {
$this->resetBeyondMinimum($sample, $beyondMinimum);
}
}
}
/**
* @param array $sample
* @param array $beyondMinimum
*/
private function resetBeyondMinimum(array &$sample, array $beyondMinimum)
{
foreach ($beyondMinimum as $index) {
$sample[$index] = 0;
}
}
/**
* @param int $samplesCount
*
* @return array
*/
private function getBeyondMinimumIndexes(int $samplesCount)
{
$indexes = [];
foreach ($this->frequencies as $token => $frequency) {
if (($frequency / $samplesCount) < $this->minDF) {
$indexes[] = $this->getTokenIndex($token);
}
}
return $indexes;
}
}
@@ -0,0 +1,197 @@
<?php
declare(strict_types=1);
namespace Phpml\Helper;
trait OneVsRest
{
/**
* @var array
*/
protected $classifiers = [];
/**
* All provided training targets' labels.
*
* @var array
*/
protected $allLabels = [];
/**
* @var array
*/
protected $costValues = [];
/**
* Train a binary classifier in the OvR style
*
* @param array $samples
* @param array $targets
*/
public function train(array $samples, array $targets)
{
// Clears previous stuff.
$this->reset();
return $this->trainBylabel($samples, $targets);
}
/**
* @param array $samples
* @param array $targets
* @param array $allLabels All training set labels
* @return void
*/
protected function trainByLabel(array $samples, array $targets, array $allLabels = [])
{
// Overwrites the current value if it exist. $allLabels must be provided for each partialTrain run.
if (!empty($allLabels)) {
$this->allLabels = $allLabels;
} else {
$this->allLabels = array_keys(array_count_values($targets));
}
sort($this->allLabels, SORT_STRING);
// If there are only two targets, then there is no need to perform OvR
if (count($this->allLabels) == 2) {
// Init classifier if required.
if (empty($this->classifiers)) {
$this->classifiers[0] = $this->getClassifierCopy();
}
$this->classifiers[0]->trainBinary($samples, $targets, $this->allLabels);
} else {
// Train a separate classifier for each label and memorize them
foreach ($this->allLabels as $label) {
// Init classifier if required.
if (empty($this->classifiers[$label])) {
$this->classifiers[$label] = $this->getClassifierCopy();
}
list($binarizedTargets, $classifierLabels) = $this->binarizeTargets($targets, $label);
$this->classifiers[$label]->trainBinary($samples, $binarizedTargets, $classifierLabels);
}
}
// If the underlying classifier is capable of giving the cost values
// during the training, then assign it to the relevant variable
// Adding just the first classifier cost values to avoid complex average calculations.
$classifierref = reset($this->classifiers);
if (method_exists($classifierref, 'getCostValues')) {
$this->costValues = $classifierref->getCostValues();
}
}
/**
* Resets the classifier and the vars internally used by OneVsRest to create multiple classifiers.
*/
public function reset()
{
$this->classifiers = [];
$this->allLabels = [];
$this->costValues = [];
$this->resetBinary();
}
/**
* Returns an instance of the current class after cleaning up OneVsRest stuff.
*
* @return \Phpml\Estimator
*/
protected function getClassifierCopy()
{
// Clone the current classifier, so that
// we don't mess up its variables while training
// multiple instances of this classifier
$classifier = clone $this;
$classifier->reset();
return $classifier;
}
/**
* Groups all targets into two groups: Targets equal to
* the given label and the others
*
* $targets is not passed by reference nor contains objects so this method
* changes will not affect the caller $targets array.
*
* @param array $targets
* @param mixed $label
* @return array Binarized targets and target's labels
*/
private function binarizeTargets($targets, $label)
{
$notLabel = "not_$label";
foreach ($targets as $key => $target) {
$targets[$key] = $target == $label ? $label : $notLabel;
}
$labels = [$label, $notLabel];
return [$targets, $labels];
}
/**
* @param array $sample
*
* @return mixed
*/
protected function predictSample(array $sample)
{
if (count($this->allLabels) == 2) {
return $this->classifiers[0]->predictSampleBinary($sample);
}
$probs = [];
foreach ($this->classifiers as $label => $predictor) {
$probs[$label] = $predictor->predictProbability($sample, $label);
}
arsort($probs, SORT_NUMERIC);
return key($probs);
}
/**
* Each classifier should implement this method instead of train(samples, targets)
*
* @param array $samples
* @param array $targets
* @param array $labels
*/
abstract protected function trainBinary(array $samples, array $targets, array $labels);
/**
* To be overwritten by OneVsRest classifiers.
*
* @return void
*/
abstract protected function resetBinary();
/**
* Each classifier that make use of OvR approach should be able to
* return a probability for a sample to belong to the given label.
*
* @param array $sample
*
* @return mixed
*/
abstract protected function predictProbability(array $sample, string $label);
/**
* Each classifier should implement this method instead of predictSample()
*
* @param array $sample
*
* @return mixed
*/
abstract protected function predictSampleBinary(array $sample);
}
@@ -0,0 +1,361 @@
<?php
declare(strict_types=1);
namespace Phpml\Helper\Optimizer;
/**
* Conjugate Gradient method to solve a non-linear f(x) with respect to unknown x
* See https://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method)
*
* The method applied below is explained in the below document in a practical manner
* - http://web.cs.iastate.edu/~cs577/handouts/conjugate-gradient.pdf
*
* However it is compliant with the general Conjugate Gradient method with
* Fletcher-Reeves update method. Note that, the f(x) is assumed to be one-dimensional
* and one gradient is utilized for all dimensions in the given data.
*/
class ConjugateGradient extends GD
{
/**
* @param array $samples
* @param array $targets
* @param \Closure $gradientCb
*
* @return array
*/
public function runOptimization(array $samples, array $targets, \Closure $gradientCb)
{
$this->samples = $samples;
$this->targets = $targets;
$this->gradientCb = $gradientCb;
$this->sampleCount = count($samples);
$this->costValues = [];
$d = mp::muls($this->gradient($this->theta), -1);
for ($i=0; $i < $this->maxIterations; $i++) {
// Obtain α that minimizes f(θ + α.d)
$alpha = $this->getAlpha(array_sum($d));
// θ(k+1) = θ(k) + α.d
$thetaNew = $this->getNewTheta($alpha, $d);
// β = ||∇f(x(k+1))||² ||∇f(x(k))||²
$beta = $this->getBeta($thetaNew);
// d(k+1) =–∇f(x(k+1)) + β(k).d(k)
$d = $this->getNewDirection($thetaNew, $beta, $d);
// Save values for the next iteration
$oldTheta = $this->theta;
$this->costValues[] = $this->cost($thetaNew);
$this->theta = $thetaNew;
if ($this->enableEarlyStop && $this->earlyStop($oldTheta)) {
break;
}
}
$this->clear();
return $this->theta;
}
/**
* Executes the callback function for the problem and returns
* sum of the gradient for all samples & targets.
*
* @param array $theta
*
* @return float
*/
protected function gradient(array $theta)
{
list($_, $gradient, $_) = parent::gradient($theta);
return $gradient;
}
/**
* Returns the value of f(x) for given solution
*
* @param array $theta
*
* @return float
*/
protected function cost(array $theta)
{
list($cost, $_, $_) = parent::gradient($theta);
return array_sum($cost) / $this->sampleCount;
}
/**
* Calculates alpha that minimizes the function f(θ + α.d)
* by performing a line search that does not rely upon the derivation.
*
* There are several alternatives for this function. For now, we
* prefer a method inspired from the bisection method for its simplicity.
* This algorithm attempts to find an optimum alpha value between 0.0001 and 0.01
*
* Algorithm as follows:
* a) Probe a small alpha (0.0001) and calculate cost function
* b) Probe a larger alpha (0.01) and calculate cost function
* b-1) If cost function decreases, continue enlarging alpha
* b-2) If cost function increases, take the midpoint and try again
*
* @param float $d
*
* @return array
*/
protected function getAlpha(float $d)
{
$small = 0.0001 * $d;
$large = 0.01 * $d;
// Obtain θ + α.d for two initial values, x0 and x1
$x0 = mp::adds($this->theta, $small);
$x1 = mp::adds($this->theta, $large);
$epsilon = 0.0001;
$iteration = 0;
do {
$fx1 = $this->cost($x1);
$fx0 = $this->cost($x0);
// If the difference between two values is small enough
// then break the loop
if (abs($fx1 - $fx0) <= $epsilon) {
break;
}
if ($fx1 < $fx0) {
$x0 = $x1;
$x1 = mp::adds($x1, 0.01); // Enlarge second
} else {
$x1 = mp::divs(mp::add($x1, $x0), 2.0);
} // Get to the midpoint
$error = $fx1 / $this->dimensions;
} while ($error <= $epsilon || $iteration++ < 10);
// Return α = θ / d
if ($d == 0) {
return $x1[0] - $this->theta[0];
}
return ($x1[0] - $this->theta[0]) / $d;
}
/**
* Calculates new set of solutions with given alpha (for each θ(k)) and
* gradient direction.
*
* θ(k+1) = θ(k) + α.d
*
* @param float $alpha
* @param array $d
*
* return array
*/
protected function getNewTheta(float $alpha, array $d)
{
$theta = $this->theta;
for ($i=0; $i < $this->dimensions + 1; $i++) {
if ($i == 0) {
$theta[$i] += $alpha * array_sum($d);
} else {
$sum = 0.0;
foreach ($this->samples as $si => $sample) {
$sum += $sample[$i - 1] * $d[$si] * $alpha;
}
$theta[$i] += $sum;
}
}
return $theta;
}
/**
* Calculates new beta (β) for given set of solutions by using
* FletcherReeves method.
*
* β = ||f(x(k+1))||² ||f(x(k))||²
*
* See:
* R. Fletcher and C. M. Reeves, "Function minimization by conjugate gradients", Comput. J. 7 (1964), 149154.
*
* @param array $newTheta
*
* @return float
*/
protected function getBeta(array $newTheta)
{
$dNew = array_sum($this->gradient($newTheta));
$dOld = array_sum($this->gradient($this->theta)) + 1e-100;
return $dNew ** 2 / $dOld ** 2;
}
/**
* Calculates the new conjugate direction
*
* d(k+1) =–∇f(x(k+1)) + β(k).d(k)
*
* @param array $theta
* @param float $beta
* @param array $d
*
* @return array
*/
protected function getNewDirection(array $theta, float $beta, array $d)
{
$grad = $this->gradient($theta);
return mp::add(mp::muls($grad, -1), mp::muls($d, $beta));
}
}
/**
* Handles element-wise vector operations between vector-vector
* and vector-scalar variables
*/
class mp
{
/**
* Element-wise <b>multiplication</b> of two vectors of the same size
*
* @param array $m1
* @param array $m2
*
* @return array
*/
public static function mul(array $m1, array $m2)
{
$res = [];
foreach ($m1 as $i => $val) {
$res[] = $val * $m2[$i];
}
return $res;
}
/**
* Element-wise <b>division</b> of two vectors of the same size
*
* @param array $m1
* @param array $m2
*
* @return array
*/
public static function div(array $m1, array $m2)
{
$res = [];
foreach ($m1 as $i => $val) {
$res[] = $val / $m2[$i];
}
return $res;
}
/**
* Element-wise <b>addition</b> of two vectors of the same size
*
* @param array $m1
* @param array $m2
*
* @return array
*/
public static function add(array $m1, array $m2, $mag = 1)
{
$res = [];
foreach ($m1 as $i => $val) {
$res[] = $val + $mag * $m2[$i];
}
return $res;
}
/**
* Element-wise <b>subtraction</b> of two vectors of the same size
*
* @param array $m1
* @param array $m2
*
* @return array
*/
public static function sub(array $m1, array $m2)
{
return self::add($m1, $m2, -1);
}
/**
* Element-wise <b>multiplication</b> of a vector with a scalar
*
* @param array $m1
* @param float $m2
*
* @return array
*/
public static function muls(array $m1, float $m2)
{
$res = [];
foreach ($m1 as $val) {
$res[] = $val * $m2;
}
return $res;
}
/**
* Element-wise <b>division</b> of a vector with a scalar
*
* @param array $m1
* @param float $m2
*
* @return array
*/
public static function divs(array $m1, float $m2)
{
$res = [];
foreach ($m1 as $val) {
$res[] = $val / ($m2 + 1e-32);
}
return $res;
}
/**
* Element-wise <b>addition</b> of a vector with a scalar
*
* @param array $m1
* @param float $m2
*
* @return array
*/
public static function adds(array $m1, float $m2, $mag = 1)
{
$res = [];
foreach ($m1 as $val) {
$res[] = $val + $mag * $m2;
}
return $res;
}
/**
* Element-wise <b>subtraction</b> of a vector with a scalar
*
* @param array $m1
* @param float $m2
*
* @return array
*/
public static function subs(array $m1, array $m2)
{
return self::adds($m1, $m2, -1);
}
}
@@ -0,0 +1,121 @@
<?php
declare(strict_types=1);
namespace Phpml\Helper\Optimizer;
/**
* Batch version of Gradient Descent to optimize the weights
* of a classifier given samples, targets and the objective function to minimize
*/
class GD extends StochasticGD
{
/**
* Number of samples given
*
* @var int
*/
protected $sampleCount = null;
/**
* @param array $samples
* @param array $targets
* @param \Closure $gradientCb
*
* @return array
*/
public function runOptimization(array $samples, array $targets, \Closure $gradientCb)
{
$this->samples = $samples;
$this->targets = $targets;
$this->gradientCb = $gradientCb;
$this->sampleCount = count($this->samples);
// Batch learning is executed:
$currIter = 0;
$this->costValues = [];
while ($this->maxIterations > $currIter++) {
$theta = $this->theta;
// Calculate update terms for each sample
list($errors, $updates, $totalPenalty) = $this->gradient($theta);
$this->updateWeightsWithUpdates($updates, $totalPenalty);
$this->costValues[] = array_sum($errors)/$this->sampleCount;
if ($this->earlyStop($theta)) {
break;
}
}
$this->clear();
return $this->theta;
}
/**
* Calculates gradient, cost function and penalty term for each sample
* then returns them as an array of values
*
* @param array $theta
*
* @return array
*/
protected function gradient(array $theta)
{
$costs = [];
$gradient= [];
$totalPenalty = 0;
foreach ($this->samples as $index => $sample) {
$target = $this->targets[$index];
$result = ($this->gradientCb)($theta, $sample, $target);
list($cost, $grad, $penalty) = array_pad($result, 3, 0);
$costs[] = $cost;
$gradient[]= $grad;
$totalPenalty += $penalty;
}
$totalPenalty /= $this->sampleCount;
return [$costs, $gradient, $totalPenalty];
}
/**
* @param array $updates
* @param float $penalty
*/
protected function updateWeightsWithUpdates(array $updates, float $penalty)
{
// Updates all weights at once
for ($i=0; $i <= $this->dimensions; $i++) {
if ($i == 0) {
$this->theta[0] -= $this->learningRate * array_sum($updates);
} else {
$col = array_column($this->samples, $i - 1);
$error = 0;
foreach ($col as $index => $val) {
$error += $val * $updates[$index];
}
$this->theta[$i] -= $this->learningRate *
($error + $penalty * $this->theta[$i]);
}
}
}
/**
* Clears the optimizer internal vars after the optimization process.
*
* @return void
*/
protected function clear()
{
$this->sampleCount = null;
parent::clear();
}
}
@@ -0,0 +1,61 @@
<?php
declare(strict_types=1);
namespace Phpml\Helper\Optimizer;
abstract class Optimizer
{
/**
* Unknown variables to be found
*
* @var array
*/
protected $theta;
/**
* Number of dimensions
*
* @var int
*/
protected $dimensions;
/**
* Inits a new instance of Optimizer for the given number of dimensions
*
* @param int $dimensions
*/
public function __construct(int $dimensions)
{
$this->dimensions = $dimensions;
// Inits the weights randomly
$this->theta = [];
for ($i=0; $i < $this->dimensions; $i++) {
$this->theta[] = rand() / (float) getrandmax();
}
}
/**
* Sets the weights manually
*
* @param array $theta
*/
public function setInitialTheta(array $theta)
{
if (count($theta) != $this->dimensions) {
throw new \Exception("Number of values in the weights array should be $this->dimensions");
}
$this->theta = $theta;
return $this;
}
/**
* Executes the optimization with the given samples & targets
* and returns the weights
*
*/
abstract protected function runOptimization(array $samples, array $targets, \Closure $gradientCb);
}
@@ -0,0 +1,285 @@
<?php
declare(strict_types=1);
namespace Phpml\Helper\Optimizer;
/**
* Stochastic Gradient Descent optimization method
* to find a solution for the equation A.ϴ = y where
* A (samples) and y (targets) are known and ϴ is unknown.
*/
class StochasticGD extends Optimizer
{
/**
* A (samples)
*
* @var array
*/
protected $samples = [];
/**
* y (targets)
*
* @var array
*/
protected $targets = [];
/**
* Callback function to get the gradient and cost value
* for a specific set of theta (ϴ) and a pair of sample & target
*
* @var \Closure
*/
protected $gradientCb = null;
/**
* Maximum number of iterations used to train the model
*
* @var int
*/
protected $maxIterations = 1000;
/**
* Learning rate is used to control the speed of the optimization.<br>
*
* Larger values of lr may overshoot the optimum or even cause divergence
* while small values slows down the convergence and increases the time
* required for the training
*
* @var float
*/
protected $learningRate = 0.001;
/**
* Minimum amount of change in the weights and error values
* between iterations that needs to be obtained to continue the training
*
* @var float
*/
protected $threshold = 1e-4;
/**
* Enable/Disable early stopping by checking the weight & cost values
* to see whether they changed large enough to continue the optimization
*
* @var bool
*/
protected $enableEarlyStop = true;
/**
* List of values obtained by evaluating the cost function at each iteration
* of the algorithm
*
* @var array
*/
protected $costValues= [];
/**
* Initializes the SGD optimizer for the given number of dimensions
*
* @param int $dimensions
*/
public function __construct(int $dimensions)
{
// Add one more dimension for the bias
parent::__construct($dimensions + 1);
$this->dimensions = $dimensions;
}
/**
* Sets minimum value for the change in the theta values
* between iterations to continue the iterations.<br>
*
* If change in the theta is less than given value then the
* algorithm will stop training
*
* @param float $threshold
*
* @return $this
*/
public function setChangeThreshold(float $threshold = 1e-5)
{
$this->threshold = $threshold;
return $this;
}
/**
* Enable/Disable early stopping by checking at each iteration
* whether changes in theta or cost value are not large enough
*
* @param bool $enable
*
* @return $this
*/
public function setEarlyStop(bool $enable = true)
{
$this->enableEarlyStop = $enable;
return $this;
}
/**
* @param float $learningRate
*
* @return $this
*/
public function setLearningRate(float $learningRate)
{
$this->learningRate = $learningRate;
return $this;
}
/**
* @param int $maxIterations
*
* @return $this
*/
public function setMaxIterations(int $maxIterations)
{
$this->maxIterations = $maxIterations;
return $this;
}
/**
* Optimization procedure finds the unknow variables for the equation A.ϴ = y
* for the given samples (A) and targets (y).<br>
*
* The cost function to minimize and the gradient of the function are to be
* handled by the callback function provided as the third parameter of the method.
*
* @param array $samples
* @param array $targets
* @param \Closure $gradientCb
*
* @return array
*/
public function runOptimization(array $samples, array $targets, \Closure $gradientCb)
{
$this->samples = $samples;
$this->targets = $targets;
$this->gradientCb = $gradientCb;
$currIter = 0;
$bestTheta = null;
$bestScore = 0.0;
$bestWeightIter = 0;
$this->costValues = [];
while ($this->maxIterations > $currIter++) {
$theta = $this->theta;
// Update the guess
$cost = $this->updateTheta();
// Save the best theta in the "pocket" so that
// any future set of theta worse than this will be disregarded
if ($bestTheta == null || $cost <= $bestScore) {
$bestTheta = $theta;
$bestScore = $cost;
$bestWeightIter = $currIter;
}
// Add the cost value for this iteration to the list
$this->costValues[] = $cost;
// Check for early stop
if ($this->enableEarlyStop && $this->earlyStop($theta)) {
break;
}
}
$this->clear();
// Solution in the pocket is better than or equal to the last state
// so, we use this solution
return $this->theta = $bestTheta;
}
/**
* @return float
*/
protected function updateTheta()
{
$jValue = 0.0;
$theta = $this->theta;
foreach ($this->samples as $index => $sample) {
$target = $this->targets[$index];
$result = ($this->gradientCb)($theta, $sample, $target);
list($error, $gradient, $penalty) = array_pad($result, 3, 0);
// Update bias
$this->theta[0] -= $this->learningRate * $gradient;
// Update other values
for ($i=1; $i <= $this->dimensions; $i++) {
$this->theta[$i] -= $this->learningRate *
($gradient * $sample[$i - 1] + $penalty * $this->theta[$i]);
}
// Sum error rate
$jValue += $error;
}
return $jValue / count($this->samples);
}
/**
* Checks if the optimization is not effective enough and can be stopped
* in case large enough changes in the solution do not happen
*
* @param array $oldTheta
*
* @return boolean
*/
protected function earlyStop($oldTheta)
{
// Check for early stop: No change larger than threshold (default 1e-5)
$diff = array_map(
function ($w1, $w2) {
return abs($w1 - $w2) > $this->threshold ? 1 : 0;
},
$oldTheta, $this->theta);
if (array_sum($diff) == 0) {
return true;
}
// Check if the last two cost values are almost the same
$costs = array_slice($this->costValues, -2);
if (count($costs) == 2 && abs($costs[1] - $costs[0]) < $this->threshold) {
return true;
}
return false;
}
/**
* Returns the list of cost values for each iteration executed in
* last run of the optimization
*
* @return array
*/
public function getCostValues()
{
return $this->costValues;
}
/**
* Clears the optimizer internal vars after the optimization process.
*
* @return void
*/
protected function clear()
{
$this->samples = [];
$this->targets = [];
$this->gradientCb = null;
}
}
@@ -0,0 +1,34 @@
<?php
declare(strict_types=1);
namespace Phpml\Helper;
trait Predictable
{
/**
* @param array $samples
*
* @return mixed
*/
public function predict(array $samples)
{
if (!is_array($samples[0])) {
$predicted = $this->predictSample($samples);
} else {
$predicted = [];
foreach ($samples as $index => $sample) {
$predicted[$index] = $this->predictSample($sample);
}
}
return $predicted;
}
/**
* @param array $sample
*
* @return mixed
*/
abstract protected function predictSample(array $sample);
}
@@ -0,0 +1,28 @@
<?php
declare(strict_types=1);
namespace Phpml\Helper;
trait Trainable
{
/**
* @var array
*/
private $samples = [];
/**
* @var array
*/
private $targets = [];
/**
* @param array $samples
* @param array $targets
*/
public function train(array $samples, array $targets)
{
$this->samples = array_merge($this->samples, $samples);
$this->targets = array_merge($this->targets, $targets);
}
}
@@ -0,0 +1,16 @@
<?php
declare(strict_types=1);
namespace Phpml;
interface IncrementalEstimator
{
/**
* @param array $samples
* @param array $targets
* @param array $labels
*/
public function partialTrain(array $samples, array $targets, array $labels = []);
}
@@ -0,0 +1,16 @@
<?php
declare(strict_types=1);
namespace Phpml\Math;
interface Distance
{
/**
* @param array $a
* @param array $b
*
* @return float
*/
public function distance(array $a, array $b): float;
}
@@ -0,0 +1,35 @@
<?php
declare(strict_types=1);
namespace Phpml\Math\Distance;
use Phpml\Exception\InvalidArgumentException;
use Phpml\Math\Distance;
class Chebyshev implements Distance
{
/**
* @param array $a
* @param array $b
*
* @return float
*
* @throws InvalidArgumentException
*/
public function distance(array $a, array $b): float
{
if (count($a) !== count($b)) {
throw InvalidArgumentException::arraySizeNotMatch();
}
$differences = [];
$count = count($a);
for ($i = 0; $i < $count; ++$i) {
$differences[] = abs($a[$i] - $b[$i]);
}
return max($differences);
}
}
@@ -0,0 +1,47 @@
<?php
declare(strict_types=1);
namespace Phpml\Math\Distance;
use Phpml\Exception\InvalidArgumentException;
use Phpml\Math\Distance;
class Euclidean implements Distance
{
/**
* @param array $a
* @param array $b
*
* @return float
*
* @throws InvalidArgumentException
*/
public function distance(array $a, array $b): float
{
if (count($a) !== count($b)) {
throw InvalidArgumentException::arraySizeNotMatch();
}
$distance = 0;
foreach ($a as $i => $val) {
$distance += ($val - $b[$i]) ** 2;
}
return sqrt((float) $distance);
}
/**
* Square of Euclidean distance
*
* @param array $a
* @param array $b
*
* @return float
*/
public function sqDistance(array $a, array $b): float
{
return $this->distance($a, $b) ** 2;
}
}
@@ -0,0 +1,35 @@
<?php
declare(strict_types=1);
namespace Phpml\Math\Distance;
use Phpml\Exception\InvalidArgumentException;
use Phpml\Math\Distance;
class Manhattan implements Distance
{
/**
* @param array $a
* @param array $b
*
* @return float
*
* @throws InvalidArgumentException
*/
public function distance(array $a, array $b): float
{
if (count($a) !== count($b)) {
throw InvalidArgumentException::arraySizeNotMatch();
}
$distance = 0;
$count = count($a);
for ($i = 0; $i < $count; ++$i) {
$distance += abs($a[$i] - $b[$i]);
}
return $distance;
}
}
@@ -0,0 +1,48 @@
<?php
declare(strict_types=1);
namespace Phpml\Math\Distance;
use Phpml\Exception\InvalidArgumentException;
use Phpml\Math\Distance;
class Minkowski implements Distance
{
/**
* @var float
*/
private $lambda;
/**
* @param float $lambda
*/
public function __construct(float $lambda = 3.0)
{
$this->lambda = $lambda;
}
/**
* @param array $a
* @param array $b
*
* @return float
*
* @throws InvalidArgumentException
*/
public function distance(array $a, array $b): float
{
if (count($a) !== count($b)) {
throw InvalidArgumentException::arraySizeNotMatch();
}
$distance = 0;
$count = count($a);
for ($i = 0; $i < $count; ++$i) {
$distance += pow(abs($a[$i] - $b[$i]), $this->lambda);
}
return (float)pow($distance, 1 / $this->lambda);
}
}
@@ -0,0 +1,16 @@
<?php
declare(strict_types=1);
namespace Phpml\Math;
interface Kernel
{
/**
* @param float $a
* @param float $b
*
* @return float
*/
public function compute($a, $b);
}
@@ -0,0 +1,39 @@
<?php
declare(strict_types=1);
namespace Phpml\Math\Kernel;
use Phpml\Math\Kernel;
use Phpml\Math\Product;
class RBF implements Kernel
{
/**
* @var float
*/
private $gamma;
/**
* @param float $gamma
*/
public function __construct(float $gamma)
{
$this->gamma = $gamma;
}
/**
* @param float $a
* @param float $b
*
* @return float
*/
public function compute($a, $b)
{
$score = 2 * Product::scalar($a, $b);
$squares = Product::scalar($a, $a) + Product::scalar($b, $b);
$result = exp(-$this->gamma * ($squares - $score));
return $result;
}
}
@@ -0,0 +1,890 @@
<?php
declare(strict_types=1);
/**
*
* Class to obtain eigenvalues and eigenvectors of a real matrix.
*
* If A is symmetric, then A = V*D*V' where the eigenvalue matrix D
* is diagonal and the eigenvector matrix V is orthogonal (i.e.
* A = V.times(D.times(V.transpose())) and V.times(V.transpose())
* equals the identity matrix).
*
* If A is not symmetric, then the eigenvalue matrix D is block diagonal
* with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
* lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
* columns of V represent the eigenvectors in the sense that A*V = V*D,
* i.e. A.times(V) equals V.times(D). The matrix V may be badly
* conditioned, or even singular, so the validity of the equation
* A = V*D*inverse(V) depends upon V.cond().
*
* @author Paul Meagher
* @license PHP v3.0
* @version 1.1
*
* Slightly changed to adapt the original code to PHP-ML library
* @date 2017/04/11
* @author Mustafa Karabulut
*/
namespace Phpml\Math\LinearAlgebra;
use Phpml\Math\Matrix;
class EigenvalueDecomposition
{
/**
* Row and column dimension (square matrix).
* @var int
*/
private $n;
/**
* Internal symmetry flag.
* @var int
*/
private $issymmetric;
/**
* Arrays for internal storage of eigenvalues.
* @var array
*/
private $d = [];
private $e = [];
/**
* Array for internal storage of eigenvectors.
* @var array
*/
private $V = [];
/**
* Array for internal storage of nonsymmetric Hessenberg form.
* @var array
*/
private $H = [];
/**
* Working storage for nonsymmetric algorithm.
* @var array
*/
private $ort;
/**
* Used for complex scalar division.
* @var float
*/
private $cdivr;
private $cdivi;
/**
* Symmetric Householder reduction to tridiagonal form.
*/
private function tred2()
{
// This is derived from the Algol procedures tred2 by
// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.
$this->d = $this->V[$this->n-1];
// Householder reduction to tridiagonal form.
for ($i = $this->n-1; $i > 0; --$i) {
$i_ = $i -1;
// Scale to avoid under/overflow.
$h = $scale = 0.0;
$scale += array_sum(array_map('abs', $this->d));
if ($scale == 0.0) {
$this->e[$i] = $this->d[$i_];
$this->d = array_slice($this->V[$i_], 0, $i_);
for ($j = 0; $j < $i; ++$j) {
$this->V[$j][$i] = $this->V[$i][$j] = 0.0;
}
} else {
// Generate Householder vector.
for ($k = 0; $k < $i; ++$k) {
$this->d[$k] /= $scale;
$h += pow($this->d[$k], 2);
}
$f = $this->d[$i_];
$g = sqrt($h);
if ($f > 0) {
$g = -$g;
}
$this->e[$i] = $scale * $g;
$h = $h - $f * $g;
$this->d[$i_] = $f - $g;
for ($j = 0; $j < $i; ++$j) {
$this->e[$j] = 0.0;
}
// Apply similarity transformation to remaining columns.
for ($j = 0; $j < $i; ++$j) {
$f = $this->d[$j];
$this->V[$j][$i] = $f;
$g = $this->e[$j] + $this->V[$j][$j] * $f;
for ($k = $j+1; $k <= $i_; ++$k) {
$g += $this->V[$k][$j] * $this->d[$k];
$this->e[$k] += $this->V[$k][$j] * $f;
}
$this->e[$j] = $g;
}
$f = 0.0;
for ($j = 0; $j < $i; ++$j) {
if ($h === 0) {
$h = 1e-20;
}
$this->e[$j] /= $h;
$f += $this->e[$j] * $this->d[$j];
}
$hh = $f / (2 * $h);
for ($j=0; $j < $i; ++$j) {
$this->e[$j] -= $hh * $this->d[$j];
}
for ($j = 0; $j < $i; ++$j) {
$f = $this->d[$j];
$g = $this->e[$j];
for ($k = $j; $k <= $i_; ++$k) {
$this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
}
$this->d[$j] = $this->V[$i-1][$j];
$this->V[$i][$j] = 0.0;
}
}
$this->d[$i] = $h;
}
// Accumulate transformations.
for ($i = 0; $i < $this->n-1; ++$i) {
$this->V[$this->n-1][$i] = $this->V[$i][$i];
$this->V[$i][$i] = 1.0;
$h = $this->d[$i+1];
if ($h != 0.0) {
for ($k = 0; $k <= $i; ++$k) {
$this->d[$k] = $this->V[$k][$i+1] / $h;
}
for ($j = 0; $j <= $i; ++$j) {
$g = 0.0;
for ($k = 0; $k <= $i; ++$k) {
$g += $this->V[$k][$i+1] * $this->V[$k][$j];
}
for ($k = 0; $k <= $i; ++$k) {
$this->V[$k][$j] -= $g * $this->d[$k];
}
}
}
for ($k = 0; $k <= $i; ++$k) {
$this->V[$k][$i+1] = 0.0;
}
}
$this->d = $this->V[$this->n-1];
$this->V[$this->n-1] = array_fill(0, $j, 0.0);
$this->V[$this->n-1][$this->n-1] = 1.0;
$this->e[0] = 0.0;
}
/**
* Symmetric tridiagonal QL algorithm.
*
* This is derived from the Algol procedures tql2, by
* Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
* Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
* Fortran subroutine in EISPACK.
*/
private function tql2()
{
for ($i = 1; $i < $this->n; ++$i) {
$this->e[$i-1] = $this->e[$i];
}
$this->e[$this->n-1] = 0.0;
$f = 0.0;
$tst1 = 0.0;
$eps = pow(2.0, -52.0);
for ($l = 0; $l < $this->n; ++$l) {
// Find small subdiagonal element
$tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
$m = $l;
while ($m < $this->n) {
if (abs($this->e[$m]) <= $eps * $tst1) {
break;
}
++$m;
}
// If m == l, $this->d[l] is an eigenvalue,
// otherwise, iterate.
if ($m > $l) {
$iter = 0;
do {
// Could check iteration count here.
$iter += 1;
// Compute implicit shift
$g = $this->d[$l];
$p = ($this->d[$l+1] - $g) / (2.0 * $this->e[$l]);
$r = hypot($p, 1.0);
if ($p < 0) {
$r *= -1;
}
$this->d[$l] = $this->e[$l] / ($p + $r);
$this->d[$l+1] = $this->e[$l] * ($p + $r);
$dl1 = $this->d[$l+1];
$h = $g - $this->d[$l];
for ($i = $l + 2; $i < $this->n; ++$i) {
$this->d[$i] -= $h;
}
$f += $h;
// Implicit QL transformation.
$p = $this->d[$m];
$c = 1.0;
$c2 = $c3 = $c;
$el1 = $this->e[$l + 1];
$s = $s2 = 0.0;
for ($i = $m-1; $i >= $l; --$i) {
$c3 = $c2;
$c2 = $c;
$s2 = $s;
$g = $c * $this->e[$i];
$h = $c * $p;
$r = hypot($p, $this->e[$i]);
$this->e[$i+1] = $s * $r;
$s = $this->e[$i] / $r;
$c = $p / $r;
$p = $c * $this->d[$i] - $s * $g;
$this->d[$i+1] = $h + $s * ($c * $g + $s * $this->d[$i]);
// Accumulate transformation.
for ($k = 0; $k < $this->n; ++$k) {
$h = $this->V[$k][$i+1];
$this->V[$k][$i+1] = $s * $this->V[$k][$i] + $c * $h;
$this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
}
}
$p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
$this->e[$l] = $s * $p;
$this->d[$l] = $c * $p;
// Check for convergence.
} while (abs($this->e[$l]) > $eps * $tst1);
}
$this->d[$l] = $this->d[$l] + $f;
$this->e[$l] = 0.0;
}
// Sort eigenvalues and corresponding vectors.
for ($i = 0; $i < $this->n - 1; ++$i) {
$k = $i;
$p = $this->d[$i];
for ($j = $i+1; $j < $this->n; ++$j) {
if ($this->d[$j] < $p) {
$k = $j;
$p = $this->d[$j];
}
}
if ($k != $i) {
$this->d[$k] = $this->d[$i];
$this->d[$i] = $p;
for ($j = 0; $j < $this->n; ++$j) {
$p = $this->V[$j][$i];
$this->V[$j][$i] = $this->V[$j][$k];
$this->V[$j][$k] = $p;
}
}
}
}
/**
* Nonsymmetric reduction to Hessenberg form.
*
* This is derived from the Algol procedures orthes and ortran,
* by Martin and Wilkinson, Handbook for Auto. Comp.,
* Vol.ii-Linear Algebra, and the corresponding
* Fortran subroutines in EISPACK.
*/
private function orthes()
{
$low = 0;
$high = $this->n-1;
for ($m = $low+1; $m <= $high-1; ++$m) {
// Scale column.
$scale = 0.0;
for ($i = $m; $i <= $high; ++$i) {
$scale = $scale + abs($this->H[$i][$m-1]);
}
if ($scale != 0.0) {
// Compute Householder transformation.
$h = 0.0;
for ($i = $high; $i >= $m; --$i) {
$this->ort[$i] = $this->H[$i][$m-1] / $scale;
$h += $this->ort[$i] * $this->ort[$i];
}
$g = sqrt($h);
if ($this->ort[$m] > 0) {
$g *= -1;
}
$h -= $this->ort[$m] * $g;
$this->ort[$m] -= $g;
// Apply Householder similarity transformation
// H = (I -u * u' / h) * H * (I -u * u') / h)
for ($j = $m; $j < $this->n; ++$j) {
$f = 0.0;
for ($i = $high; $i >= $m; --$i) {
$f += $this->ort[$i] * $this->H[$i][$j];
}
$f /= $h;
for ($i = $m; $i <= $high; ++$i) {
$this->H[$i][$j] -= $f * $this->ort[$i];
}
}
for ($i = 0; $i <= $high; ++$i) {
$f = 0.0;
for ($j = $high; $j >= $m; --$j) {
$f += $this->ort[$j] * $this->H[$i][$j];
}
$f = $f / $h;
for ($j = $m; $j <= $high; ++$j) {
$this->H[$i][$j] -= $f * $this->ort[$j];
}
}
$this->ort[$m] = $scale * $this->ort[$m];
$this->H[$m][$m-1] = $scale * $g;
}
}
// Accumulate transformations (Algol's ortran).
for ($i = 0; $i < $this->n; ++$i) {
for ($j = 0; $j < $this->n; ++$j) {
$this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
}
}
for ($m = $high-1; $m >= $low+1; --$m) {
if ($this->H[$m][$m-1] != 0.0) {
for ($i = $m+1; $i <= $high; ++$i) {
$this->ort[$i] = $this->H[$i][$m-1];
}
for ($j = $m; $j <= $high; ++$j) {
$g = 0.0;
for ($i = $m; $i <= $high; ++$i) {
$g += $this->ort[$i] * $this->V[$i][$j];
}
// Double division avoids possible underflow
$g = ($g / $this->ort[$m]) / $this->H[$m][$m-1];
for ($i = $m; $i <= $high; ++$i) {
$this->V[$i][$j] += $g * $this->ort[$i];
}
}
}
}
}
/**
* Performs complex division.
*/
private function cdiv($xr, $xi, $yr, $yi)
{
if (abs($yr) > abs($yi)) {
$r = $yi / $yr;
$d = $yr + $r * $yi;
$this->cdivr = ($xr + $r * $xi) / $d;
$this->cdivi = ($xi - $r * $xr) / $d;
} else {
$r = $yr / $yi;
$d = $yi + $r * $yr;
$this->cdivr = ($r * $xr + $xi) / $d;
$this->cdivi = ($r * $xi - $xr) / $d;
}
}
/**
* Nonsymmetric reduction from Hessenberg to real Schur form.
*
* Code is derived from the Algol procedure hqr2,
* by Martin and Wilkinson, Handbook for Auto. Comp.,
* Vol.ii-Linear Algebra, and the corresponding
* Fortran subroutine in EISPACK.
*/
private function hqr2()
{
// Initialize
$nn = $this->n;
$n = $nn - 1;
$low = 0;
$high = $nn - 1;
$eps = pow(2.0, -52.0);
$exshift = 0.0;
$p = $q = $r = $s = $z = 0;
// Store roots isolated by balanc and compute matrix norm
$norm = 0.0;
for ($i = 0; $i < $nn; ++$i) {
if (($i < $low) or ($i > $high)) {
$this->d[$i] = $this->H[$i][$i];
$this->e[$i] = 0.0;
}
for ($j = max($i-1, 0); $j < $nn; ++$j) {
$norm = $norm + abs($this->H[$i][$j]);
}
}
// Outer loop over eigenvalue index
$iter = 0;
while ($n >= $low) {
// Look for single small sub-diagonal element
$l = $n;
while ($l > $low) {
$s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]);
if ($s == 0.0) {
$s = $norm;
}
if (abs($this->H[$l][$l-1]) < $eps * $s) {
break;
}
--$l;
}
// Check for convergence
// One root found
if ($l == $n) {
$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
$this->d[$n] = $this->H[$n][$n];
$this->e[$n] = 0.0;
--$n;
$iter = 0;
// Two roots found
} elseif ($l == $n-1) {
$w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
$p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0;
$q = $p * $p + $w;
$z = sqrt(abs($q));
$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
$this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift;
$x = $this->H[$n][$n];
// Real pair
if ($q >= 0) {
if ($p >= 0) {
$z = $p + $z;
} else {
$z = $p - $z;
}
$this->d[$n-1] = $x + $z;
$this->d[$n] = $this->d[$n-1];
if ($z != 0.0) {
$this->d[$n] = $x - $w / $z;
}
$this->e[$n-1] = 0.0;
$this->e[$n] = 0.0;
$x = $this->H[$n][$n-1];
$s = abs($x) + abs($z);
$p = $x / $s;
$q = $z / $s;
$r = sqrt($p * $p + $q * $q);
$p = $p / $r;
$q = $q / $r;
// Row modification
for ($j = $n-1; $j < $nn; ++$j) {
$z = $this->H[$n-1][$j];
$this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j];
$this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
}
// Column modification
for ($i = 0; $i <= $n; ++$i) {
$z = $this->H[$i][$n-1];
$this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n];
$this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
}
// Accumulate transformations
for ($i = $low; $i <= $high; ++$i) {
$z = $this->V[$i][$n-1];
$this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n];
$this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
}
// Complex pair
} else {
$this->d[$n-1] = $x + $p;
$this->d[$n] = $x + $p;
$this->e[$n-1] = $z;
$this->e[$n] = -$z;
}
$n = $n - 2;
$iter = 0;
// No convergence yet
} else {
// Form shift
$x = $this->H[$n][$n];
$y = 0.0;
$w = 0.0;
if ($l < $n) {
$y = $this->H[$n-1][$n-1];
$w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
}
// Wilkinson's original ad hoc shift
if ($iter == 10) {
$exshift += $x;
for ($i = $low; $i <= $n; ++$i) {
$this->H[$i][$i] -= $x;
}
$s = abs($this->H[$n][$n-1]) + abs($this->H[$n-1][$n-2]);
$x = $y = 0.75 * $s;
$w = -0.4375 * $s * $s;
}
// MATLAB's new ad hoc shift
if ($iter == 30) {
$s = ($y - $x) / 2.0;
$s = $s * $s + $w;
if ($s > 0) {
$s = sqrt($s);
if ($y < $x) {
$s = -$s;
}
$s = $x - $w / (($y - $x) / 2.0 + $s);
for ($i = $low; $i <= $n; ++$i) {
$this->H[$i][$i] -= $s;
}
$exshift += $s;
$x = $y = $w = 0.964;
}
}
// Could check iteration count here.
$iter = $iter + 1;
// Look for two consecutive small sub-diagonal elements
$m = $n - 2;
while ($m >= $l) {
$z = $this->H[$m][$m];
$r = $x - $z;
$s = $y - $z;
$p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1];
$q = $this->H[$m+1][$m+1] - $z - $r - $s;
$r = $this->H[$m+2][$m+1];
$s = abs($p) + abs($q) + abs($r);
$p = $p / $s;
$q = $q / $s;
$r = $r / $s;
if ($m == $l) {
break;
}
if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) <
$eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) {
break;
}
--$m;
}
for ($i = $m + 2; $i <= $n; ++$i) {
$this->H[$i][$i-2] = 0.0;
if ($i > $m+2) {
$this->H[$i][$i-3] = 0.0;
}
}
// Double QR step involving rows l:n and columns m:n
for ($k = $m; $k <= $n-1; ++$k) {
$notlast = ($k != $n-1);
if ($k != $m) {
$p = $this->H[$k][$k-1];
$q = $this->H[$k+1][$k-1];
$r = ($notlast ? $this->H[$k+2][$k-1] : 0.0);
$x = abs($p) + abs($q) + abs($r);
if ($x != 0.0) {
$p = $p / $x;
$q = $q / $x;
$r = $r / $x;
}
}
if ($x == 0.0) {
break;
}
$s = sqrt($p * $p + $q * $q + $r * $r);
if ($p < 0) {
$s = -$s;
}
if ($s != 0) {
if ($k != $m) {
$this->H[$k][$k-1] = -$s * $x;
} elseif ($l != $m) {
$this->H[$k][$k-1] = -$this->H[$k][$k-1];
}
$p = $p + $s;
$x = $p / $s;
$y = $q / $s;
$z = $r / $s;
$q = $q / $p;
$r = $r / $p;
// Row modification
for ($j = $k; $j < $nn; ++$j) {
$p = $this->H[$k][$j] + $q * $this->H[$k+1][$j];
if ($notlast) {
$p = $p + $r * $this->H[$k+2][$j];
$this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z;
}
$this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
$this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y;
}
// Column modification
for ($i = 0; $i <= min($n, $k+3); ++$i) {
$p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1];
if ($notlast) {
$p = $p + $z * $this->H[$i][$k+2];
$this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r;
}
$this->H[$i][$k] = $this->H[$i][$k] - $p;
$this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q;
}
// Accumulate transformations
for ($i = $low; $i <= $high; ++$i) {
$p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1];
if ($notlast) {
$p = $p + $z * $this->V[$i][$k+2];
$this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r;
}
$this->V[$i][$k] = $this->V[$i][$k] - $p;
$this->V[$i][$k+1] = $this->V[$i][$k+1] - $p * $q;
}
} // ($s != 0)
} // k loop
} // check convergence
} // while ($n >= $low)
// Backsubstitute to find vectors of upper triangular form
if ($norm == 0.0) {
return;
}
for ($n = $nn-1; $n >= 0; --$n) {
$p = $this->d[$n];
$q = $this->e[$n];
// Real vector
if ($q == 0) {
$l = $n;
$this->H[$n][$n] = 1.0;
for ($i = $n-1; $i >= 0; --$i) {
$w = $this->H[$i][$i] - $p;
$r = 0.0;
for ($j = $l; $j <= $n; ++$j) {
$r = $r + $this->H[$i][$j] * $this->H[$j][$n];
}
if ($this->e[$i] < 0.0) {
$z = $w;
$s = $r;
} else {
$l = $i;
if ($this->e[$i] == 0.0) {
if ($w != 0.0) {
$this->H[$i][$n] = -$r / $w;
} else {
$this->H[$i][$n] = -$r / ($eps * $norm);
}
// Solve real equations
} else {
$x = $this->H[$i][$i+1];
$y = $this->H[$i+1][$i];
$q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
$t = ($x * $s - $z * $r) / $q;
$this->H[$i][$n] = $t;
if (abs($x) > abs($z)) {
$this->H[$i+1][$n] = (-$r - $w * $t) / $x;
} else {
$this->H[$i+1][$n] = (-$s - $y * $t) / $z;
}
}
// Overflow control
$t = abs($this->H[$i][$n]);
if (($eps * $t) * $t > 1) {
for ($j = $i; $j <= $n; ++$j) {
$this->H[$j][$n] = $this->H[$j][$n] / $t;
}
}
}
}
// Complex vector
} elseif ($q < 0) {
$l = $n-1;
// Last vector component imaginary so matrix is triangular
if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) {
$this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1];
$this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1];
} else {
$this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q);
$this->H[$n-1][$n-1] = $this->cdivr;
$this->H[$n-1][$n] = $this->cdivi;
}
$this->H[$n][$n-1] = 0.0;
$this->H[$n][$n] = 1.0;
for ($i = $n-2; $i >= 0; --$i) {
// double ra,sa,vr,vi;
$ra = 0.0;
$sa = 0.0;
for ($j = $l; $j <= $n; ++$j) {
$ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1];
$sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
}
$w = $this->H[$i][$i] - $p;
if ($this->e[$i] < 0.0) {
$z = $w;
$r = $ra;
$s = $sa;
} else {
$l = $i;
if ($this->e[$i] == 0) {
$this->cdiv(-$ra, -$sa, $w, $q);
$this->H[$i][$n-1] = $this->cdivr;
$this->H[$i][$n] = $this->cdivi;
} else {
// Solve complex equations
$x = $this->H[$i][$i+1];
$y = $this->H[$i+1][$i];
$vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
$vi = ($this->d[$i] - $p) * 2.0 * $q;
if ($vr == 0.0 & $vi == 0.0) {
$vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
}
$this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
$this->H[$i][$n-1] = $this->cdivr;
$this->H[$i][$n] = $this->cdivi;
if (abs($x) > (abs($z) + abs($q))) {
$this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x;
$this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x;
} else {
$this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q);
$this->H[$i+1][$n-1] = $this->cdivr;
$this->H[$i+1][$n] = $this->cdivi;
}
}
// Overflow control
$t = max(abs($this->H[$i][$n-1]), abs($this->H[$i][$n]));
if (($eps * $t) * $t > 1) {
for ($j = $i; $j <= $n; ++$j) {
$this->H[$j][$n-1] = $this->H[$j][$n-1] / $t;
$this->H[$j][$n] = $this->H[$j][$n] / $t;
}
}
} // end else
} // end for
} // end else for complex case
} // end for
// Vectors of isolated roots
for ($i = 0; $i < $nn; ++$i) {
if ($i < $low | $i > $high) {
for ($j = $i; $j < $nn; ++$j) {
$this->V[$i][$j] = $this->H[$i][$j];
}
}
}
// Back transformation to get eigenvectors of original matrix
for ($j = $nn-1; $j >= $low; --$j) {
for ($i = $low; $i <= $high; ++$i) {
$z = 0.0;
for ($k = $low; $k <= min($j, $high); ++$k) {
$z = $z + $this->V[$i][$k] * $this->H[$k][$j];
}
$this->V[$i][$j] = $z;
}
}
} // end hqr2
/**
* Constructor: Check for symmetry, then construct the eigenvalue decomposition
*
* @param array $Arg
*/
public function __construct(array $Arg)
{
$this->A = $Arg;
$this->n = count($Arg[0]);
$issymmetric = true;
for ($j = 0; ($j < $this->n) & $issymmetric; ++$j) {
for ($i = 0; ($i < $this->n) & $issymmetric; ++$i) {
$issymmetric = ($this->A[$i][$j] == $this->A[$j][$i]);
}
}
if ($issymmetric) {
$this->V = $this->A;
// Tridiagonalize.
$this->tred2();
// Diagonalize.
$this->tql2();
} else {
$this->H = $this->A;
$this->ort = [];
// Reduce to Hessenberg form.
$this->orthes();
// Reduce Hessenberg to real Schur form.
$this->hqr2();
}
}
/**
* Return the eigenvector matrix
*
* @access public
* @return array
*/
public function getEigenvectors()
{
$vectors = $this->V;
// Always return the eigenvectors of length 1.0
$vectors = new Matrix($vectors);
$vectors = array_map(function ($vect) {
$sum = 0;
for ($i=0; $i<count($vect); $i++) {
$sum += $vect[$i] ** 2;
}
$sum = sqrt($sum);
for ($i=0; $i<count($vect); $i++) {
$vect[$i] /= $sum;
}
return $vect;
}, $vectors->transpose()->toArray());
return $vectors;
}
/**
* Return the real parts of the eigenvalues<br>
* d = real(diag(D));
*
* @return array
*/
public function getRealEigenvalues()
{
return $this->d;
}
/**
* Return the imaginary parts of the eigenvalues <br>
* d = imag(diag(D))
*
* @return array
*/
public function getImagEigenvalues()
{
return $this->e;
}
/**
* Return the block diagonal eigenvalue matrix
*
* @return array
*/
public function getDiagonalEigenvalues()
{
for ($i = 0; $i < $this->n; ++$i) {
$D[$i] = array_fill(0, $this->n, 0.0);
$D[$i][$i] = $this->d[$i];
if ($this->e[$i] == 0) {
continue;
}
$o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
$D[$i][$o] = $this->e[$i];
}
return $D;
}
} // class EigenvalueDecomposition
@@ -0,0 +1,385 @@
<?php
declare(strict_types=1);
namespace Phpml\Math;
use Phpml\Exception\InvalidArgumentException;
use Phpml\Exception\MatrixException;
class Matrix
{
/**
* @var array
*/
private $matrix;
/**
* @var int
*/
private $rows;
/**
* @var int
*/
private $columns;
/**
* @var float
*/
private $determinant;
/**
* @param array $matrix
* @param bool $validate
*
* @throws InvalidArgumentException
*/
public function __construct(array $matrix, bool $validate = true)
{
// When a row vector is given
if (!is_array($matrix[0])) {
$this->rows = 1;
$this->columns = count($matrix);
$matrix = [$matrix];
} else {
$this->rows = count($matrix);
$this->columns = count($matrix[0]);
}
if ($validate) {
for ($i = 0; $i < $this->rows; ++$i) {
if (count($matrix[$i]) !== $this->columns) {
throw InvalidArgumentException::matrixDimensionsDidNotMatch();
}
}
}
$this->matrix = $matrix;
}
/**
* @param array $array
*
* @return Matrix
*/
public static function fromFlatArray(array $array)
{
$matrix = [];
foreach ($array as $value) {
$matrix[] = [$value];
}
return new self($matrix);
}
/**
* @return array
*/
public function toArray()
{
return $this->matrix;
}
/**
* @return float
*/
public function toScalar()
{
return $this->matrix[0][0];
}
/**
* @return int
*/
public function getRows()
{
return $this->rows;
}
/**
* @return int
*/
public function getColumns()
{
return $this->columns;
}
/**
* @param $column
*
* @return array
*
* @throws MatrixException
*/
public function getColumnValues($column)
{
if ($column >= $this->columns) {
throw MatrixException::columnOutOfRange();
}
return array_column($this->matrix, $column);
}
/**
* @return float|int
*
* @throws MatrixException
*/
public function getDeterminant()
{
if ($this->determinant) {
return $this->determinant;
}
if (!$this->isSquare()) {
throw MatrixException::notSquareMatrix();
}
return $this->determinant = $this->calculateDeterminant();
}
/**
* @return float|int
*
* @throws MatrixException
*/
private function calculateDeterminant()
{
$determinant = 0;
if ($this->rows == 1 && $this->columns == 1) {
$determinant = $this->matrix[0][0];
} elseif ($this->rows == 2 && $this->columns == 2) {
$determinant =
$this->matrix[0][0] * $this->matrix[1][1] -
$this->matrix[0][1] * $this->matrix[1][0];
} else {
for ($j = 0; $j < $this->columns; ++$j) {
$subMatrix = $this->crossOut(0, $j);
$minor = $this->matrix[0][$j] * $subMatrix->getDeterminant();
$determinant += fmod((float) $j, 2.0) == 0 ? $minor : -$minor;
}
}
return $determinant;
}
/**
* @return bool
*/
public function isSquare()
{
return $this->columns === $this->rows;
}
/**
* @return Matrix
*/
public function transpose()
{
if ($this->rows == 1) {
$matrix = array_map(function ($el) {
return [$el];
}, $this->matrix[0]);
} else {
$matrix = array_map(null, ...$this->matrix);
}
return new self($matrix, false);
}
/**
* @param Matrix $matrix
*
* @return Matrix
*
* @throws InvalidArgumentException
*/
public function multiply(Matrix $matrix)
{
if ($this->columns != $matrix->getRows()) {
throw InvalidArgumentException::inconsistentMatrixSupplied();
}
$product = [];
$multiplier = $matrix->toArray();
for ($i = 0; $i < $this->rows; ++$i) {
$columns = $matrix->getColumns();
for ($j = 0; $j < $columns; ++$j) {
$product[$i][$j] = 0;
for ($k = 0; $k < $this->columns; ++$k) {
$product[$i][$j] += $this->matrix[$i][$k] * $multiplier[$k][$j];
}
}
}
return new self($product, false);
}
/**
* @param $value
*
* @return Matrix
*/
public function divideByScalar($value)
{
$newMatrix = [];
for ($i = 0; $i < $this->rows; ++$i) {
for ($j = 0; $j < $this->columns; ++$j) {
$newMatrix[$i][$j] = $this->matrix[$i][$j] / $value;
}
}
return new self($newMatrix, false);
}
/**
* @param $value
*
* @return Matrix
*/
public function multiplyByScalar($value)
{
$newMatrix = [];
for ($i = 0; $i < $this->rows; ++$i) {
for ($j = 0; $j < $this->columns; ++$j) {
$newMatrix[$i][$j] = $this->matrix[$i][$j] * $value;
}
}
return new self($newMatrix, false);
}
/**
* Element-wise addition of the matrix with another one
*
* @param Matrix $other
*/
public function add(Matrix $other)
{
return $this->_add($other);
}
/**
* Element-wise subtracting of another matrix from this one
*
* @param Matrix $other
*/
public function subtract(Matrix $other)
{
return $this->_add($other, -1);
}
/**
* Element-wise addition or substraction depending on the given sign parameter
*
* @param Matrix $other
* @param type $sign
*/
protected function _add(Matrix $other, $sign = 1)
{
$a1 = $this->toArray();
$a2 = $other->toArray();
$newMatrix = [];
for ($i=0; $i < $this->rows; $i++) {
for ($k=0; $k < $this->columns; $k++) {
$newMatrix[$i][$k] = $a1[$i][$k] + $sign * $a2[$i][$k];
}
}
return new Matrix($newMatrix, false);
}
/**
* @return Matrix
*
* @throws MatrixException
*/
public function inverse()
{
if (!$this->isSquare()) {
throw MatrixException::notSquareMatrix();
}
if ($this->isSingular()) {
throw MatrixException::singularMatrix();
}
$newMatrix = [];
for ($i = 0; $i < $this->rows; ++$i) {
for ($j = 0; $j < $this->columns; ++$j) {
$minor = $this->crossOut($i, $j)->getDeterminant();
$newMatrix[$i][$j] = fmod((float) ($i + $j), 2.0) == 0 ? $minor : -$minor;
}
}
$cofactorMatrix = new self($newMatrix, false);
return $cofactorMatrix->transpose()->divideByScalar($this->getDeterminant());
}
/**
* @param int $row
* @param int $column
*
* @return Matrix
*/
public function crossOut(int $row, int $column)
{
$newMatrix = [];
$r = 0;
for ($i = 0; $i < $this->rows; ++$i) {
$c = 0;
if ($row != $i) {
for ($j = 0; $j < $this->columns; ++$j) {
if ($column != $j) {
$newMatrix[$r][$c] = $this->matrix[$i][$j];
++$c;
}
}
++$r;
}
}
return new self($newMatrix, false);
}
/**
* @return bool
*/
public function isSingular() : bool
{
return 0 == $this->getDeterminant();
}
/**
* Returns the transpose of given array
*
* @param array $array
*
* @return array
*/
public static function transposeArray(array $array)
{
return (new Matrix($array, false))->transpose()->toArray();
}
/**
* Returns the dot product of two arrays<br>
* Matrix::dot(x, y) ==> x.y'
*
* @param array $array1
* @param array $array2
*
* @return array
*/
public static function dot(array $array1, array $array2)
{
$m1 = new Matrix($array1, false);
$m2 = new Matrix($array2, false);
return $m1->multiply($m2->transpose())->toArray()[0];
}
}
@@ -0,0 +1,26 @@
<?php
declare(strict_types=1);
namespace Phpml\Math;
class Product
{
/**
* @param array $a
* @param array $b
*
* @return mixed
*/
public static function scalar(array $a, array $b)
{
$product = 0;
foreach ($a as $index => $value) {
if (is_numeric($value) && is_numeric($b[$index])) {
$product += $value * $b[$index];
}
}
return $product;
}
}
@@ -0,0 +1,211 @@
<?php
declare(strict_types=1);
namespace Phpml\Math;
class Set implements \IteratorAggregate
{
/**
* @var string[]|int[]|float[]
*/
private $elements;
/**
* @param string[]|int[]|float[] $elements
*/
public function __construct(array $elements = [])
{
$this->elements = self::sanitize($elements);
}
/**
* Creates the union of A and B.
*
* @param Set $a
* @param Set $b
*
* @return Set
*/
public static function union(Set $a, Set $b) : Set
{
return new self(array_merge($a->toArray(), $b->toArray()));
}
/**
* Creates the intersection of A and B.
*
* @param Set $a
* @param Set $b
*
* @return Set
*/
public static function intersection(Set $a, Set $b) : Set
{
return new self(array_intersect($a->toArray(), $b->toArray()));
}
/**
* Creates the difference of A and B.
*
* @param Set $a
* @param Set $b
*
* @return Set
*/
public static function difference(Set $a, Set $b) : Set
{
return new self(array_diff($a->toArray(), $b->toArray()));
}
/**
* Creates the Cartesian product of A and B.
*
* @param Set $a
* @param Set $b
*
* @return Set[]
*/
public static function cartesian(Set $a, Set $b) : array
{
$cartesian = [];
foreach ($a as $multiplier) {
foreach ($b as $multiplicand) {
$cartesian[] = new self(array_merge([$multiplicand], [$multiplier]));
}
}
return $cartesian;
}
/**
* Creates the power set of A.
*
* @param Set $a
*
* @return Set[]
*/
public static function power(Set $a) : array
{
$power = [new self()];
foreach ($a as $multiplicand) {
foreach ($power as $multiplier) {
$power[] = new self(array_merge([$multiplicand], $multiplier->toArray()));
}
}
return $power;
}
/**
* Removes duplicates and rewrites index.
*
* @param string[]|int[]|float[] $elements
*
* @return string[]|int[]|float[]
*/
private static function sanitize(array $elements) : array
{
sort($elements, SORT_ASC);
return array_values(array_unique($elements, SORT_ASC));
}
/**
* @param string|int|float $element
*
* @return Set
*/
public function add($element) : Set
{
return $this->addAll([$element]);
}
/**
* @param string[]|int[]|float[] $elements
*
* @return Set
*/
public function addAll(array $elements) : Set
{
$this->elements = self::sanitize(array_merge($this->elements, $elements));
return $this;
}
/**
* @param string|int|float $element
*
* @return Set
*/
public function remove($element) : Set
{
return $this->removeAll([$element]);
}
/**
* @param string[]|int[]|float[] $elements
*
* @return Set
*/
public function removeAll(array $elements) : Set
{
$this->elements = self::sanitize(array_diff($this->elements, $elements));
return $this;
}
/**
* @param string|int|float $element
*
* @return bool
*/
public function contains($element) : bool
{
return $this->containsAll([$element]);
}
/**
* @param string[]|int[]|float[] $elements
*
* @return bool
*/
public function containsAll(array $elements) : bool
{
return !array_diff($elements, $this->elements);
}
/**
* @return string[]|int[]|float[]
*/
public function toArray() : array
{
return $this->elements;
}
/**
* @return \ArrayIterator
*/
public function getIterator() : \ArrayIterator
{
return new \ArrayIterator($this->elements);
}
/**
* @return bool
*/
public function isEmpty() : bool
{
return $this->cardinality() == 0;
}
/**
* @return int
*/
public function cardinality() : int
{
return count($this->elements);
}
}
@@ -0,0 +1,45 @@
<?php
declare(strict_types=1);
namespace Phpml\Math\Statistic;
use Phpml\Exception\InvalidArgumentException;
class Correlation
{
/**
* @param array|int[]|float[] $x
* @param array|int[]|float[] $y
*
* @return float
*
* @throws InvalidArgumentException
*/
public static function pearson(array $x, array $y)
{
if (count($x) !== count($y)) {
throw InvalidArgumentException::arraySizeNotMatch();
}
$count = count($x);
$meanX = Mean::arithmetic($x);
$meanY = Mean::arithmetic($y);
$axb = 0;
$a2 = 0;
$b2 = 0;
for ($i = 0; $i < $count; ++$i) {
$a = $x[$i] - $meanX;
$b = $y[$i] - $meanY;
$axb += ($a * $b);
$a2 += pow($a, 2);
$b2 += pow($b, 2);
}
$corr = $axb / sqrt((float) ($a2 * $b2));
return $corr;
}
}
@@ -0,0 +1,155 @@
<?php
declare(strict_types=1);
namespace Phpml\Math\Statistic;
use Phpml\Exception\InvalidArgumentException;
class Covariance
{
/**
* Calculates covariance from two given arrays, x and y, respectively
*
* @param array $x
* @param array $y
* @param bool $sample
* @param float $meanX
* @param float $meanY
*
* @return float
*
* @throws InvalidArgumentException
*/
public static function fromXYArrays(array $x, array $y, $sample = true, float $meanX = null, float $meanY = null)
{
if (empty($x) || empty($y)) {
throw InvalidArgumentException::arrayCantBeEmpty();
}
$n = count($x);
if ($sample && $n === 1) {
throw InvalidArgumentException::arraySizeToSmall(2);
}
if ($meanX === null) {
$meanX = Mean::arithmetic($x);
}
if ($meanY === null) {
$meanY = Mean::arithmetic($y);
}
$sum = 0.0;
foreach ($x as $index => $xi) {
$yi = $y[$index];
$sum += ($xi - $meanX) * ($yi - $meanY);
}
if ($sample) {
--$n;
}
return $sum / $n;
}
/**
* Calculates covariance of two dimensions, i and k in the given data.
*
* @param array $data
* @param int $i
* @param int $k
* @param type $sample
* @param int $n
* @param float $meanX
* @param float $meanY
*/
public static function fromDataset(array $data, int $i, int $k, $sample = true, float $meanX = null, float $meanY = null)
{
if (empty($data)) {
throw InvalidArgumentException::arrayCantBeEmpty();
}
$n = count($data);
if ($sample && $n === 1) {
throw InvalidArgumentException::arraySizeToSmall(2);
}
if ($i < 0 || $k < 0 || $i >= $n || $k >= $n) {
throw new \Exception("Given indices i and k do not match with the dimensionality of data");
}
if ($meanX === null || $meanY === null) {
$x = array_column($data, $i);
$y = array_column($data, $k);
$meanX = Mean::arithmetic($x);
$meanY = Mean::arithmetic($y);
$sum = 0.0;
foreach ($x as $index => $xi) {
$yi = $y[$index];
$sum += ($xi - $meanX) * ($yi - $meanY);
}
} else {
// In the case, whole dataset given along with dimension indices, i and k,
// we would like to avoid getting column data with array_column and operate
// over this extra copy of column data for memory efficiency purposes.
//
// Instead we traverse through the whole data and get what we actually need
// without copying the data. This way, memory use will be reduced
// with a slight cost of CPU utilization.
$sum = 0.0;
foreach ($data as $row) {
$val = [];
foreach ($row as $index => $col) {
if ($index == $i) {
$val[0] = $col - $meanX;
}
if ($index == $k) {
$val[1] = $col - $meanY;
}
}
$sum += $val[0] * $val[1];
}
}
if ($sample) {
--$n;
}
return $sum / $n;
}
/**
* Returns the covariance matrix of n-dimensional data
*
* @param array $data
*
* @return array
*/
public static function covarianceMatrix(array $data, array $means = null)
{
$n = count($data[0]);
if ($means === null) {
$means = [];
for ($i=0; $i < $n; $i++) {
$means[] = Mean::arithmetic(array_column($data, $i));
}
}
$cov = [];
for ($i=0; $i < $n; $i++) {
for ($k=0; $k < $n; $k++) {
if ($i > $k) {
$cov[$i][$k] = $cov[$k][$i];
} else {
$cov[$i][$k] = Covariance::fromDataset(
$data, $i, $k, true, $means[$i], $means[$k]);
}
}
}
return $cov;
}
}
@@ -0,0 +1,60 @@
<?php
declare(strict_types=1);
namespace Phpml\Math\Statistic;
class Gaussian
{
/**
* @var float
*/
protected $mean;
/**
* @var float
*/
protected $std;
/**
* @param float $mean
* @param float $std
*/
public function __construct(float $mean, float $std)
{
$this->mean = $mean;
$this->std = $std;
}
/**
* Returns probability density of the given <i>$value</i>
*
* @param float $value
*
* @return type
*/
public function pdf(float $value)
{
// Calculate the probability density by use of normal/Gaussian distribution
// Ref: https://en.wikipedia.org/wiki/Normal_distribution
$std2 = $this->std ** 2;
$mean = $this->mean;
return exp(- (($value - $mean) ** 2) / (2 * $std2)) / sqrt(2 * $std2 * pi());
}
/**
* Returns probability density value of the given <i>$value</i> based on
* given standard deviation and the mean
*
* @param float $mean
* @param float $std
* @param float $value
*
* @return float
*/
public static function distributionPdf(float $mean, float $std, float $value)
{
$normal = new self($mean, $std);
return $normal->pdf($value);
}
}
@@ -0,0 +1,75 @@
<?php
declare(strict_types=1);
namespace Phpml\Math\Statistic;
use Phpml\Exception\InvalidArgumentException;
class Mean
{
/**
* @param array $numbers
*
* @return float
*
* @throws InvalidArgumentException
*/
public static function arithmetic(array $numbers)
{
self::checkArrayLength($numbers);
return array_sum($numbers) / count($numbers);
}
/**
* @param array $numbers
*
* @return float|mixed
*
* @throws InvalidArgumentException
*/
public static function median(array $numbers)
{
self::checkArrayLength($numbers);
$count = count($numbers);
$middleIndex = (int)floor($count / 2);
sort($numbers, SORT_NUMERIC);
$median = $numbers[$middleIndex];
if (0 === $count % 2) {
$median = ($median + $numbers[$middleIndex - 1]) / 2;
}
return $median;
}
/**
* @param array $numbers
*
* @return mixed
*
* @throws InvalidArgumentException
*/
public static function mode(array $numbers)
{
self::checkArrayLength($numbers);
$values = array_count_values($numbers);
return array_search(max($values), $values);
}
/**
* @param array $array
*
* @throws InvalidArgumentException
*/
private static function checkArrayLength(array $array)
{
if (0 == count($array)) {
throw InvalidArgumentException::arrayCantBeEmpty();
}
}
}
@@ -0,0 +1,44 @@
<?php
declare(strict_types=1);
namespace Phpml\Math\Statistic;
use Phpml\Exception\InvalidArgumentException;
class StandardDeviation
{
/**
* @param array|float[] $a
* @param bool $sample
*
* @return float
*
* @throws InvalidArgumentException
*/
public static function population(array $a, $sample = true)
{
if (empty($a)) {
throw InvalidArgumentException::arrayCantBeEmpty();
}
$n = count($a);
if ($sample && $n === 1) {
throw InvalidArgumentException::arraySizeToSmall(2);
}
$mean = Mean::arithmetic($a);
$carry = 0.0;
foreach ($a as $val) {
$d = $val - $mean;
$carry += $d * $d;
}
if ($sample) {
--$n;
}
return sqrt((float) ($carry / $n));
}
}
@@ -0,0 +1,39 @@
<?php
declare(strict_types=1);
namespace Phpml\Metric;
use Phpml\Exception\InvalidArgumentException;
class Accuracy
{
/**
* @param array $actualLabels
* @param array $predictedLabels
* @param bool $normalize
*
* @return float|int
*
* @throws InvalidArgumentException
*/
public static function score(array $actualLabels, array $predictedLabels, bool $normalize = true)
{
if (count($actualLabels) != count($predictedLabels)) {
throw InvalidArgumentException::arraySizeNotMatch();
}
$score = 0;
foreach ($actualLabels as $index => $label) {
if ($label == $predictedLabels[$index]) {
++$score;
}
}
if ($normalize) {
$score /= count($actualLabels);
}
return $score;
}
}
@@ -0,0 +1,183 @@
<?php
declare(strict_types=1);
namespace Phpml\Metric;
class ClassificationReport
{
/**
* @var array
*/
private $precision = [];
/**
* @var array
*/
private $recall = [];
/**
* @var array
*/
private $f1score = [];
/**
* @var array
*/
private $support = [];
/**
* @var array
*/
private $average = [];
/**
* @param array $actualLabels
* @param array $predictedLabels
*/
public function __construct(array $actualLabels, array $predictedLabels)
{
$truePositive = $falsePositive = $falseNegative = $this->support = self::getLabelIndexedArray($actualLabels, $predictedLabels);
foreach ($actualLabels as $index => $actual) {
$predicted = $predictedLabels[$index];
++$this->support[$actual];
if ($actual === $predicted) {
++$truePositive[$actual];
} else {
++$falsePositive[$predicted];
++$falseNegative[$actual];
}
}
$this->computeMetrics($truePositive, $falsePositive, $falseNegative);
$this->computeAverage();
}
/**
* @return array
*/
public function getPrecision()
{
return $this->precision;
}
/**
* @return array
*/
public function getRecall()
{
return $this->recall;
}
/**
* @return array
*/
public function getF1score()
{
return $this->f1score;
}
/**
* @return array
*/
public function getSupport()
{
return $this->support;
}
/**
* @return array
*/
public function getAverage()
{
return $this->average;
}
/**
* @param array $truePositive
* @param array $falsePositive
* @param array $falseNegative
*/
private function computeMetrics(array $truePositive, array $falsePositive, array $falseNegative)
{
foreach ($truePositive as $label => $tp) {
$this->precision[$label] = $this->computePrecision($tp, $falsePositive[$label]);
$this->recall[$label] = $this->computeRecall($tp, $falseNegative[$label]);
$this->f1score[$label] = $this->computeF1Score((float) $this->precision[$label], (float) $this->recall[$label]);
}
}
private function computeAverage()
{
foreach (['precision', 'recall', 'f1score'] as $metric) {
$values = array_filter($this->$metric);
if (0 == count($values)) {
$this->average[$metric] = 0.0;
continue;
}
$this->average[$metric] = array_sum($values) / count($values);
}
}
/**
* @param int $truePositive
* @param int $falsePositive
*
* @return float|string
*/
private function computePrecision(int $truePositive, int $falsePositive)
{
if (0 == ($divider = $truePositive + $falsePositive)) {
return 0.0;
}
return $truePositive / $divider;
}
/**
* @param int $truePositive
* @param int $falseNegative
*
* @return float|string
*/
private function computeRecall(int $truePositive, int $falseNegative)
{
if (0 == ($divider = $truePositive + $falseNegative)) {
return 0.0;
}
return $truePositive / $divider;
}
/**
* @param float $precision
* @param float $recall
*
* @return float
*/
private function computeF1Score(float $precision, float $recall): float
{
if (0 == ($divider = $precision + $recall)) {
return 0.0;
}
return 2.0 * (($precision * $recall) / $divider);
}
/**
* @param array $actualLabels
* @param array $predictedLabels
*
* @return array
*/
private static function getLabelIndexedArray(array $actualLabels, array $predictedLabels): array
{
$labels = array_values(array_unique(array_merge($actualLabels, $predictedLabels)));
sort($labels);
$labels = array_combine($labels, array_fill(0, count($labels), 0));
return $labels;
}
}
@@ -0,0 +1,71 @@
<?php
declare(strict_types=1);
namespace Phpml\Metric;
class ConfusionMatrix
{
/**
* @param array $actualLabels
* @param array $predictedLabels
* @param array $labels
*
* @return array
*/
public static function compute(array $actualLabels, array $predictedLabels, array $labels = null): array
{
$labels = $labels ? array_flip($labels) : self::getUniqueLabels($actualLabels);
$matrix = self::generateMatrixWithZeros($labels);
foreach ($actualLabels as $index => $actual) {
$predicted = $predictedLabels[$index];
if (!isset($labels[$actual]) || !isset($labels[$predicted])) {
continue;
}
if ($predicted === $actual) {
$row = $column = $labels[$actual];
} else {
$row = $labels[$actual];
$column = $labels[$predicted];
}
$matrix[$row][$column] += 1;
}
return $matrix;
}
/**
* @param array $labels
*
* @return array
*/
private static function generateMatrixWithZeros(array $labels): array
{
$count = count($labels);
$matrix = [];
for ($i = 0; $i < $count; ++$i) {
$matrix[$i] = array_fill(0, $count, 0);
}
return $matrix;
}
/**
* @param array $labels
*
* @return array
*/
private static function getUniqueLabels(array $labels): array
{
$labels = array_values(array_unique($labels));
sort($labels);
$labels = array_flip($labels);
return $labels;
}
}
@@ -0,0 +1,54 @@
<?php
declare(strict_types=1);
namespace Phpml;
use Phpml\Exception\SerializeException;
use Phpml\Exception\FileException;
class ModelManager
{
/**
* @param Estimator $estimator
* @param string $filepath
* @throws FileException
* @throws SerializeException
*/
public function saveToFile(Estimator $estimator, string $filepath)
{
if (!is_writable(dirname($filepath))) {
throw FileException::cantSaveFile(basename($filepath));
}
$serialized = serialize($estimator);
if (empty($serialized)) {
throw SerializeException::cantSerialize(get_type($estimator));
}
$result = file_put_contents($filepath, $serialized, LOCK_EX);
if ($result === false) {
throw FileException::cantSaveFile(basename($filepath));
}
}
/**
* @param string $filepath
* @return Estimator
* @throws FileException
* @throws SerializeException
*/
public function restoreFromFile(string $filepath) : Estimator
{
if (!file_exists($filepath) || !is_readable($filepath)) {
throw FileException::cantOpenFile(basename($filepath));
}
$object = unserialize(file_get_contents($filepath));
if ($object === false) {
throw SerializeException::cantUnserialize(basename($filepath));
}
return $object;
}
}
@@ -0,0 +1,15 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork;
interface ActivationFunction
{
/**
* @param float|int $value
*
* @return float
*/
public function compute($value): float;
}
@@ -0,0 +1,20 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\ActivationFunction;
use Phpml\NeuralNetwork\ActivationFunction;
class BinaryStep implements ActivationFunction
{
/**
* @param float|int $value
*
* @return float
*/
public function compute($value): float
{
return $value >= 0 ? 1.0 : 0.0;
}
}
@@ -0,0 +1,20 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\ActivationFunction;
use Phpml\NeuralNetwork\ActivationFunction;
class Gaussian implements ActivationFunction
{
/**
* @param float|int $value
*
* @return float
*/
public function compute($value): float
{
return exp(-pow($value, 2));
}
}
@@ -0,0 +1,33 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\ActivationFunction;
use Phpml\NeuralNetwork\ActivationFunction;
class HyperbolicTangent implements ActivationFunction
{
/**
* @var float
*/
private $beta;
/**
* @param float $beta
*/
public function __construct($beta = 1.0)
{
$this->beta = $beta;
}
/**
* @param float|int $value
*
* @return float
*/
public function compute($value): float
{
return tanh($this->beta * $value);
}
}
@@ -0,0 +1,33 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\ActivationFunction;
use Phpml\NeuralNetwork\ActivationFunction;
class Sigmoid implements ActivationFunction
{
/**
* @var float
*/
private $beta;
/**
* @param float $beta
*/
public function __construct($beta = 1.0)
{
$this->beta = $beta;
}
/**
* @param float|int $value
*
* @return float
*/
public function compute($value): float
{
return 1 / (1 + exp(-$this->beta * $value));
}
}
@@ -0,0 +1,65 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork;
use Phpml\Exception\InvalidArgumentException;
use Phpml\NeuralNetwork\Node\Neuron;
class Layer
{
/**
* @var Node[]
*/
private $nodes = [];
/**
* @param int $nodesNumber
* @param string $nodeClass
* @param ActivationFunction|null $activationFunction
*
* @throws InvalidArgumentException
*/
public function __construct(int $nodesNumber = 0, string $nodeClass = Neuron::class, ActivationFunction $activationFunction = null)
{
if (!in_array(Node::class, class_implements($nodeClass))) {
throw InvalidArgumentException::invalidLayerNodeClass();
}
for ($i = 0; $i < $nodesNumber; ++$i) {
$this->nodes[] = $this->createNode($nodeClass, $activationFunction);
}
}
/**
* @param string $nodeClass
* @param ActivationFunction|null $activationFunction
*
* @return Neuron
*/
private function createNode(string $nodeClass, ActivationFunction $activationFunction = null)
{
if (Neuron::class == $nodeClass) {
return new Neuron($activationFunction);
}
return new $nodeClass();
}
/**
* @param Node $node
*/
public function addNode(Node $node)
{
$this->nodes[] = $node;
}
/**
* @return Node[]
*/
public function getNodes()
{
return $this->nodes;
}
}
@@ -0,0 +1,30 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork;
interface Network
{
/**
* @param mixed $input
*
* @return self
*/
public function setInput($input);
/**
* @return array
*/
public function getOutput(): array;
/**
* @param Layer $layer
*/
public function addLayer(Layer $layer);
/**
* @return Layer[]
*/
public function getLayers(): array;
}
@@ -0,0 +1,81 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\Network;
use Phpml\NeuralNetwork\Layer;
use Phpml\NeuralNetwork\Network;
use Phpml\NeuralNetwork\Node\Input;
use Phpml\NeuralNetwork\Node\Neuron;
abstract class LayeredNetwork implements Network
{
/**
* @var Layer[]
*/
protected $layers;
/**
* @param Layer $layer
*/
public function addLayer(Layer $layer)
{
$this->layers[] = $layer;
}
/**
* @return Layer[]
*/
public function getLayers(): array
{
return $this->layers;
}
/**
* @return Layer
*/
public function getOutputLayer(): Layer
{
return $this->layers[count($this->layers) - 1];
}
/**
* @return array
*/
public function getOutput(): array
{
$result = [];
foreach ($this->getOutputLayer()->getNodes() as $neuron) {
$result[] = $neuron->getOutput();
}
return $result;
}
/**
* @param mixed $input
*
* @return $this
*/
public function setInput($input)
{
$firstLayer = $this->layers[0];
foreach ($firstLayer->getNodes() as $key => $neuron) {
if ($neuron instanceof Input) {
$neuron->setInput($input[$key]);
}
}
foreach ($this->getLayers() as $layer) {
foreach ($layer->getNodes() as $node) {
if ($node instanceof Neuron) {
$node->refresh();
}
}
}
return $this;
}
}
@@ -0,0 +1,95 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\Network;
use Phpml\Exception\InvalidArgumentException;
use Phpml\NeuralNetwork\ActivationFunction;
use Phpml\NeuralNetwork\Layer;
use Phpml\NeuralNetwork\Node\Bias;
use Phpml\NeuralNetwork\Node\Input;
use Phpml\NeuralNetwork\Node\Neuron;
use Phpml\NeuralNetwork\Node\Neuron\Synapse;
class MultilayerPerceptron extends LayeredNetwork
{
/**
* @param array $layers
* @param ActivationFunction|null $activationFunction
*
* @throws InvalidArgumentException
*/
public function __construct(array $layers, ActivationFunction $activationFunction = null)
{
if (count($layers) < 2) {
throw InvalidArgumentException::invalidLayersNumber();
}
$this->addInputLayer(array_shift($layers));
$this->addNeuronLayers($layers, $activationFunction);
$this->addBiasNodes();
$this->generateSynapses();
}
/**
* @param int $nodes
*/
private function addInputLayer(int $nodes)
{
$this->addLayer(new Layer($nodes, Input::class));
}
/**
* @param array $layers
* @param ActivationFunction|null $activationFunction
*/
private function addNeuronLayers(array $layers, ActivationFunction $activationFunction = null)
{
foreach ($layers as $neurons) {
$this->addLayer(new Layer($neurons, Neuron::class, $activationFunction));
}
}
private function generateSynapses()
{
$layersNumber = count($this->layers) - 1;
for ($i = 0; $i < $layersNumber; ++$i) {
$currentLayer = $this->layers[$i];
$nextLayer = $this->layers[$i + 1];
$this->generateLayerSynapses($nextLayer, $currentLayer);
}
}
private function addBiasNodes()
{
$biasLayers = count($this->layers) - 1;
for ($i = 0; $i < $biasLayers; ++$i) {
$this->layers[$i]->addNode(new Bias());
}
}
/**
* @param Layer $nextLayer
* @param Layer $currentLayer
*/
private function generateLayerSynapses(Layer $nextLayer, Layer $currentLayer)
{
foreach ($nextLayer->getNodes() as $nextNeuron) {
if ($nextNeuron instanceof Neuron) {
$this->generateNeuronSynapses($currentLayer, $nextNeuron);
}
}
}
/**
* @param Layer $currentLayer
* @param Neuron $nextNeuron
*/
private function generateNeuronSynapses(Layer $currentLayer, Neuron $nextNeuron)
{
foreach ($currentLayer->getNodes() as $currentNeuron) {
$nextNeuron->addSynapse(new Synapse($currentNeuron));
}
}
}
@@ -0,0 +1,13 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork;
interface Node
{
/**
* @return float
*/
public function getOutput(): float;
}
@@ -0,0 +1,18 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\Node;
use Phpml\NeuralNetwork\Node;
class Bias implements Node
{
/**
* @return float
*/
public function getOutput(): float
{
return 1.0;
}
}
@@ -0,0 +1,39 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\Node;
use Phpml\NeuralNetwork\Node;
class Input implements Node
{
/**
* @var float
*/
private $input;
/**
* @param float $input
*/
public function __construct(float $input = 0.0)
{
$this->input = $input;
}
/**
* @return float
*/
public function getOutput(): float
{
return $this->input;
}
/**
* @param float $input
*/
public function setInput(float $input)
{
$this->input = $input;
}
}
@@ -0,0 +1,75 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\Node;
use Phpml\NeuralNetwork\ActivationFunction;
use Phpml\NeuralNetwork\Node\Neuron\Synapse;
use Phpml\NeuralNetwork\Node;
class Neuron implements Node
{
/**
* @var Synapse[]
*/
protected $synapses;
/**
* @var ActivationFunction
*/
protected $activationFunction;
/**
* @var float
*/
protected $output;
/**
* @param ActivationFunction|null $activationFunction
*/
public function __construct(ActivationFunction $activationFunction = null)
{
$this->activationFunction = $activationFunction ?: new ActivationFunction\Sigmoid();
$this->synapses = [];
$this->output = 0;
}
/**
* @param Synapse $synapse
*/
public function addSynapse(Synapse $synapse)
{
$this->synapses[] = $synapse;
}
/**
* @return Synapse[]
*/
public function getSynapses()
{
return $this->synapses;
}
/**
* @return float
*/
public function getOutput(): float
{
if (0 === $this->output) {
$sum = 0;
foreach ($this->synapses as $synapse) {
$sum += $synapse->getOutput();
}
$this->output = $this->activationFunction->compute($sum);
}
return $this->output;
}
public function refresh()
{
$this->output = 0;
}
}
@@ -0,0 +1,70 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\Node\Neuron;
use Phpml\NeuralNetwork\Node;
class Synapse
{
/**
* @var float
*/
protected $weight;
/**
* @var Node
*/
protected $node;
/**
* @param Node $node
* @param float|null $weight
*/
public function __construct(Node $node, float $weight = null)
{
$this->node = $node;
$this->weight = $weight ?: $this->generateRandomWeight();
}
/**
* @return float
*/
protected function generateRandomWeight(): float
{
return 1 / random_int(5, 25) * (random_int(0, 1) ? -1 : 1);
}
/**
* @return float
*/
public function getOutput(): float
{
return $this->weight * $this->node->getOutput();
}
/**
* @param float $delta
*/
public function changeWeight($delta)
{
$this->weight += $delta;
}
/**
* @return float
*/
public function getWeight()
{
return $this->weight;
}
/**
* @return Node
*/
public function getNode()
{
return $this->node;
}
}
@@ -0,0 +1,16 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork;
interface Training
{
/**
* @param array $samples
* @param array $targets
* @param float $desiredError
* @param int $maxIterations
*/
public function train(array $samples, array $targets, float $desiredError = 0.001, int $maxIterations = 10000);
}
@@ -0,0 +1,160 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\Training;
use Phpml\NeuralNetwork\Network;
use Phpml\NeuralNetwork\Node\Neuron;
use Phpml\NeuralNetwork\Training;
use Phpml\NeuralNetwork\Training\Backpropagation\Sigma;
class Backpropagation implements Training
{
/**
* @var Network
*/
private $network;
/**
* @var int
*/
private $theta;
/**
* @var array
*/
private $sigmas;
/**
* @param Network $network
* @param int $theta
*/
public function __construct(Network $network, int $theta = 1)
{
$this->network = $network;
$this->theta = $theta;
}
/**
* @param array $samples
* @param array $targets
* @param float $desiredError
* @param int $maxIterations
*/
public function train(array $samples, array $targets, float $desiredError = 0.001, int $maxIterations = 10000)
{
for ($i = 0; $i < $maxIterations; ++$i) {
$resultsWithinError = $this->trainSamples($samples, $targets, $desiredError);
if ($resultsWithinError == count($samples)) {
break;
}
}
}
/**
* @param array $samples
* @param array $targets
* @param float $desiredError
*
* @return int
*/
private function trainSamples(array $samples, array $targets, float $desiredError): int
{
$resultsWithinError = 0;
foreach ($targets as $key => $target) {
$result = $this->network->setInput($samples[$key])->getOutput();
if ($this->isResultWithinError($result, $target, $desiredError)) {
++$resultsWithinError;
} else {
$this->trainSample($samples[$key], $target);
}
}
return $resultsWithinError;
}
/**
* @param array $sample
* @param array $target
*/
private function trainSample(array $sample, array $target)
{
$this->network->setInput($sample)->getOutput();
$this->sigmas = [];
$layers = $this->network->getLayers();
$layersNumber = count($layers);
for ($i = $layersNumber; $i > 1; --$i) {
foreach ($layers[$i - 1]->getNodes() as $key => $neuron) {
if ($neuron instanceof Neuron) {
$sigma = $this->getSigma($neuron, $target, $key, $i == $layersNumber);
foreach ($neuron->getSynapses() as $synapse) {
$synapse->changeWeight($this->theta * $sigma * $synapse->getNode()->getOutput());
}
}
}
}
}
/**
* @param Neuron $neuron
* @param array $target
* @param int $key
* @param bool $lastLayer
*
* @return float
*/
private function getSigma(Neuron $neuron, array $target, int $key, bool $lastLayer): float
{
$neuronOutput = $neuron->getOutput();
$sigma = $neuronOutput * (1 - $neuronOutput);
if ($lastLayer) {
$sigma *= ($target[$key] - $neuronOutput);
} else {
$sigma *= $this->getPrevSigma($neuron);
}
$this->sigmas[] = new Sigma($neuron, $sigma);
return $sigma;
}
/**
* @param Neuron $neuron
*
* @return float
*/
private function getPrevSigma(Neuron $neuron): float
{
$sigma = 0.0;
foreach ($this->sigmas as $neuronSigma) {
$sigma += $neuronSigma->getSigmaForNeuron($neuron);
}
return $sigma;
}
/**
* @param array $result
* @param array $target
* @param float $desiredError
*
* @return bool
*/
private function isResultWithinError(array $result, array $target, float $desiredError)
{
foreach ($target as $key => $value) {
if ($result[$key] > $value + $desiredError || $result[$key] < $value - $desiredError) {
return false;
}
}
return true;
}
}
@@ -0,0 +1,64 @@
<?php
declare(strict_types=1);
namespace Phpml\NeuralNetwork\Training\Backpropagation;
use Phpml\NeuralNetwork\Node\Neuron;
class Sigma
{
/**
* @var Neuron
*/
private $neuron;
/**
* @var float
*/
private $sigma;
/**
* @param Neuron $neuron
* @param float $sigma
*/
public function __construct(Neuron $neuron, $sigma)
{
$this->neuron = $neuron;
$this->sigma = $sigma;
}
/**
* @return Neuron
*/
public function getNeuron()
{
return $this->neuron;
}
/**
* @return float
*/
public function getSigma()
{
return $this->sigma;
}
/**
* @param Neuron $neuron
*
* @return float
*/
public function getSigmaForNeuron(Neuron $neuron): float
{
$sigma = 0.0;
foreach ($this->neuron->getSynapses() as $synapse) {
if ($synapse->getNode() == $neuron) {
$sigma += $synapse->getWeight() * $this->getSigma();
}
}
return $sigma;
}
}
@@ -0,0 +1,106 @@
<?php
declare(strict_types=1);
namespace Phpml;
class Pipeline implements Estimator
{
/**
* @var array|Transformer[]
*/
private $transformers;
/**
* @var Estimator
*/
private $estimator;
/**
* @param array|Transformer[] $transformers
* @param Estimator $estimator
*/
public function __construct(array $transformers, Estimator $estimator)
{
foreach ($transformers as $transformer) {
$this->addTransformer($transformer);
}
$this->estimator = $estimator;
}
/**
* @param Transformer $transformer
*/
public function addTransformer(Transformer $transformer)
{
$this->transformers[] = $transformer;
}
/**
* @param Estimator $estimator
*/
public function setEstimator(Estimator $estimator)
{
$this->estimator = $estimator;
}
/**
* @return array|Transformer[]
*/
public function getTransformers()
{
return $this->transformers;
}
/**
* @return Estimator
*/
public function getEstimator()
{
return $this->estimator;
}
/**
* @param array $samples
* @param array $targets
*/
public function train(array $samples, array $targets)
{
$this->fitTransformers($samples);
$this->transformSamples($samples);
$this->estimator->train($samples, $targets);
}
/**
* @param array $samples
*
* @return mixed
*/
public function predict(array $samples)
{
$this->transformSamples($samples);
return $this->estimator->predict($samples);
}
/**
* @param array $samples
*/
private function fitTransformers(array &$samples)
{
foreach ($this->transformers as $transformer) {
$transformer->fit($samples);
}
}
/**
* @param array $samples
*/
private function transformSamples(array &$samples)
{
foreach ($this->transformers as $transformer) {
$transformer->transform($samples);
}
}
}
@@ -0,0 +1,99 @@
<?php
declare(strict_types=1);
namespace Phpml\Preprocessing;
use Phpml\Preprocessing\Imputer\Strategy;
class Imputer implements Preprocessor
{
const AXIS_COLUMN = 0;
const AXIS_ROW = 1;
/**
* @var mixed
*/
private $missingValue;
/**
* @var Strategy
*/
private $strategy;
/**
* @var int
*/
private $axis;
/**
* @var
*/
private $samples;
/**
* @param mixed $missingValue
* @param Strategy $strategy
* @param int $axis
* @param array|null $samples
*/
public function __construct($missingValue, Strategy $strategy, int $axis = self::AXIS_COLUMN, array $samples = [])
{
$this->missingValue = $missingValue;
$this->strategy = $strategy;
$this->axis = $axis;
$this->samples = $samples;
}
/**
* @param array $samples
*/
public function fit(array $samples)
{
$this->samples = $samples;
}
/**
* @param array $samples
*/
public function transform(array &$samples)
{
foreach ($samples as &$sample) {
$this->preprocessSample($sample);
}
}
/**
* @param array $sample
*/
private function preprocessSample(array &$sample)
{
foreach ($sample as $column => &$value) {
if ($value === $this->missingValue) {
$value = $this->strategy->replaceValue($this->getAxis($column, $sample));
}
}
}
/**
* @param int $column
* @param array $currentSample
*
* @return array
*/
private function getAxis(int $column, array $currentSample): array
{
if (self::AXIS_ROW === $this->axis) {
return array_diff($currentSample, [$this->missingValue]);
}
$axis = [];
foreach ($this->samples as $sample) {
if ($sample[$column] !== $this->missingValue) {
$axis[] = $sample[$column];
}
}
return $axis;
}
}
@@ -0,0 +1,15 @@
<?php
declare(strict_types=1);
namespace Phpml\Preprocessing\Imputer;
interface Strategy
{
/**
* @param array $currentAxis
*
* @return mixed
*/
public function replaceValue(array $currentAxis);
}
@@ -0,0 +1,21 @@
<?php
declare(strict_types=1);
namespace Phpml\Preprocessing\Imputer\Strategy;
use Phpml\Preprocessing\Imputer\Strategy;
use Phpml\Math\Statistic\Mean;
class MeanStrategy implements Strategy
{
/**
* @param array $currentAxis
*
* @return float
*/
public function replaceValue(array $currentAxis)
{
return Mean::arithmetic($currentAxis);
}
}
@@ -0,0 +1,21 @@
<?php
declare(strict_types=1);
namespace Phpml\Preprocessing\Imputer\Strategy;
use Phpml\Preprocessing\Imputer\Strategy;
use Phpml\Math\Statistic\Mean;
class MedianStrategy implements Strategy
{
/**
* @param array $currentAxis
*
* @return float
*/
public function replaceValue(array $currentAxis)
{
return Mean::median($currentAxis);
}
}
@@ -0,0 +1,21 @@
<?php
declare(strict_types=1);
namespace Phpml\Preprocessing\Imputer\Strategy;
use Phpml\Preprocessing\Imputer\Strategy;
use Phpml\Math\Statistic\Mean;
class MostFrequentStrategy implements Strategy
{
/**
* @param array $currentAxis
*
* @return float|mixed
*/
public function replaceValue(array $currentAxis)
{
return Mean::mode($currentAxis);
}
}
@@ -0,0 +1,145 @@
<?php
declare(strict_types=1);
namespace Phpml\Preprocessing;
use Phpml\Exception\NormalizerException;
use Phpml\Math\Statistic\StandardDeviation;
use Phpml\Math\Statistic\Mean;
class Normalizer implements Preprocessor
{
const NORM_L1 = 1;
const NORM_L2 = 2;
const NORM_STD= 3;
/**
* @var int
*/
private $norm;
/**
* @var bool
*/
private $fitted = false;
/**
* @var array
*/
private $std;
/**
* @var array
*/
private $mean;
/**
* @param int $norm
*
* @throws NormalizerException
*/
public function __construct(int $norm = self::NORM_L2)
{
if (!in_array($norm, [self::NORM_L1, self::NORM_L2, self::NORM_STD])) {
throw NormalizerException::unknownNorm();
}
$this->norm = $norm;
}
/**
* @param array $samples
*/
public function fit(array $samples)
{
if ($this->fitted) {
return;
}
if ($this->norm == self::NORM_STD) {
$features = range(0, count($samples[0]) - 1);
foreach ($features as $i) {
$values = array_column($samples, $i);
$this->std[$i] = StandardDeviation::population($values);
$this->mean[$i] = Mean::arithmetic($values);
}
}
$this->fitted = true;
}
/**
* @param array $samples
*/
public function transform(array &$samples)
{
$methods = [
self::NORM_L1 => 'normalizeL1',
self::NORM_L2 => 'normalizeL2',
self::NORM_STD=> 'normalizeSTD'
];
$method = $methods[$this->norm];
$this->fit($samples);
foreach ($samples as &$sample) {
$this->$method($sample);
}
}
/**
* @param array $sample
*/
private function normalizeL1(array &$sample)
{
$norm1 = 0;
foreach ($sample as $feature) {
$norm1 += abs($feature);
}
if (0 == $norm1) {
$count = count($sample);
$sample = array_fill(0, $count, 1.0 / $count);
} else {
foreach ($sample as &$feature) {
$feature /= $norm1;
}
}
}
/**
* @param array $sample
*/
private function normalizeL2(array &$sample)
{
$norm2 = 0;
foreach ($sample as $feature) {
$norm2 += $feature * $feature;
}
$norm2 = sqrt((float)$norm2);
if (0 == $norm2) {
$sample = array_fill(0, count($sample), 1);
} else {
foreach ($sample as &$feature) {
$feature /= $norm2;
}
}
}
/**
* @param array $sample
*/
private function normalizeSTD(array &$sample)
{
foreach ($sample as $i => $val) {
if ($this->std[$i] != 0) {
$sample[$i] = ($sample[$i] - $this->mean[$i]) / $this->std[$i];
} else {
// Same value for all samples.
$sample[$i] = 0;
}
}
}
}
@@ -0,0 +1,11 @@
<?php
declare(strict_types=1);
namespace Phpml\Preprocessing;
use Phpml\Transformer;
interface Preprocessor extends Transformer
{
}

Some files were not shown because too many files have changed in this diff Show More