using MathNet.Numerics.LinearAlgebra;

namespace QuantConnect.Algorithm.CSharp
{
    /// <summary>
    /// Basic template algorithm simply initializes the date range and cash
    /// </summary>
    public class PortfolioOptimizationNumericsAlgorithm : QCAlgorithm
    {
        private const double _targetReturn = 0.1;
        private const double _riskFreeRate = 0.01;
        private double _lagrangeMultiplier;
        private double _portfolioRisk;
        private Dictionary<string, double> _mean;
        private Dictionary<string, double> _stddev;
        private Dictionary<string, double> _weights;
        private Matrix<double> R;
        private Matrix<double> S;
        private Matrix<double> Sigma;
        private Vector<double> _discountMeanVector;

        /// <summary>
        /// Initialise the data and resolution required, as well as the cash and start-end dates for your algorithm. All algorithms must initialized.
        /// </summary>
        public override void Initialize()
        {
            SetStartDate(2013, 10, 07);  //Set Start Date
            SetEndDate(2013, 10, 11);    //Set End Date
            SetCash(100000);             //Set Strategy Cash
            // Find more symbols here: http://quantconnect.com/data
            AddEquity("SPY", Resolution.Second);

            _mean = new Dictionary<string, double>
            {
                {"A", 0.04 },
                {"B", 0.08 },
                {"C", 0.12 },
                {"D", 0.15 },
            };

            _stddev = new Dictionary<string, double>
            {
                {"A", 0.07 },
                {"B", 0.12 },
                {"C", 0.18 },
                {"D", 0.26 },
            };

            _weights = _mean.ToDictionary(k => k.Key, v => 0.0);

            S = Matrix<double>.Build.DenseOfDiagonalArray(_stddev.Values.ToArray());

            R = Matrix<double>.Build.DenseOfColumnMajor(4, 4, new[] 
            {
                1.0, 0.2, 0.5, 0.3,
                0.2, 1.0, 0.7, 0.4,
                0.5, 0.7, 1.0, 0.9,
                0.3, 0.4, 0.9, 1.0
            });

            Sigma = S * R * S;

            ComputeLagrangeMultiplier();
            ComputeWeights();
            ComputePortfolioRisk();

            Log(string.Format("Lagrange Multiplier: {0,7:F4}", _lagrangeMultiplier));
            Log(string.Format("Portfolio Risk:      {0,7:P2} ", _portfolioRisk));
            foreach (var symbol in _mean.Keys)
            {
                Log(string.Format("{0}: {1,10:P2}\t{2,7:P2}\t{3,7:P2}", symbol, _weights[symbol], _mean[symbol], _stddev[symbol]));
            }
        }

        /// <summary>
        /// OnData event is the primary entry point for your algorithm. Each new data point will be pumped in here.
        /// </summary>
        /// <param name="data">Slice object keyed by symbol containing the stock data</param>
        public override void OnData(Slice data)
        {
            if (!Portfolio.Invested)
            {
                SetHoldings("SPY", 1);
                Debug("Purchased Stock");
            }
        }

        private void ComputeLagrangeMultiplier()
        {
            _discountMeanVector =
                Vector<double>.Build.DenseOfArray(_mean.Values.ToArray()) -
                Vector<double>.Build.Dense(_mean.Count, _riskFreeRate);

            var denominatorMatrix = _discountMeanVector * Sigma.Inverse() * _discountMeanVector.ToColumnMatrix();

            _lagrangeMultiplier = (_targetReturn - _riskFreeRate) / denominatorMatrix.ToArray().First();
        }

        private void ComputeWeights()
        {
            var weights = _lagrangeMultiplier * Sigma.Inverse() * _discountMeanVector.ToColumnMatrix();

            for (var i = 0; i < weights.RowCount; i++)
            {
                var kvp = _weights.ElementAt(i);
                _weights[kvp.Key] = weights.ToArray()[i, 0];
            }
        }

        private void ComputePortfolioRisk()
        {
            var weights = Vector<double>.Build.DenseOfArray(_weights.Values.ToArray());
            var portfolioVarianceMatrix = weights * Sigma * weights.ToColumnMatrix();
            _portfolioRisk = Math.Sqrt(portfolioVarianceMatrix.ToArray().First());
        }
    }
}