source: branches/3022-FastFunctionExtraction/FFX/FastFunctionExtraction.cs @ 17738

using System;
using System.Threading;
using System.Linq;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HEAL.Attic;
using HeuristicLab.Algorithms.DataAnalysis.Glmnet;
using HeuristicLab.Problems.DataAnalysis;
using System.Collections.Generic;
using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;

namespace HeuristicLab.Algorithms.DataAnalysis.FastFunctionExtraction {

    [Item(Name = "FastFunctionExtraction", Description = "An FFX algorithm.")]
    [Creatable(Category = CreatableAttribute.Categories.Algorithms, Priority = 999)]
    [StorableType("689280F7-E371-44A2-98A5-FCEDF22CA343")] // for persistence (storing your algorithm to a files or transfer to HeuristicLab.Hive
    public sealed class FastFunctionExtraction : FixedDataAnalysisAlgorithm<RegressionProblem> {

        #region constants
        private static readonly HashSet<double> exponents = new HashSet<double> { -1.0, -0.5, +0.5, +1.0 };
        private static readonly HashSet<NonlinearOperator> nonlinFuncs = new HashSet<NonlinearOperator> { NonlinearOperator.Abs, NonlinearOperator.Log, NonlinearOperator.None };
        private static readonly double minHingeThr = 0.2;
        private static readonly double maxHingeThr = 0.8;
        private static readonly int numHingeThrs = 5;

        private const string ConsiderInteractionsParameterName = "Consider Interactions";
        private const string ConsiderDenominationParameterName = "Consider Denomination";
        private const string ConsiderExponentiationParameterName = "Consider Exponentiation";
        private const string ConsiderHingeFuncsParameterName = "Consider Hinge Functions";
        private const string ConsiderNonlinearFuncsParameterName = "Consider Nonlinear Functions";
        private const string LambdaParameterName = "Elastic Net Lambda";
        private const string PenaltyParameterName = "Elastic Net Penalty";
        private const string MaxNumBasisFuncsParameterName = "Maximum Number of Basis Functions";

        #endregion

        #region parameters
        public IValueParameter<BoolValue> ConsiderInteractionsParameter {
            get { return (IValueParameter<BoolValue>)Parameters[ConsiderInteractionsParameterName]; }
        }
        public IValueParameter<BoolValue> ConsiderDenominationsParameter {
            get { return (IValueParameter<BoolValue>)Parameters[ConsiderDenominationParameterName]; }
        }
        public IValueParameter<BoolValue> ConsiderExponentiationsParameter {
            get { return (IValueParameter<BoolValue>)Parameters[ConsiderExponentiationParameterName]; }
        }
        public IValueParameter<BoolValue> ConsiderNonlinearFuncsParameter {
            get { return (IValueParameter<BoolValue>)Parameters[ConsiderNonlinearFuncsParameterName]; }
        }
        public IValueParameter<BoolValue> ConsiderHingeFuncsParameter {
            get { return (IValueParameter<BoolValue>)Parameters[ConsiderHingeFuncsParameterName]; }
        }
        public IValueParameter<DoubleValue> PenaltyParameter {
            get { return (IValueParameter<DoubleValue>)Parameters[PenaltyParameterName]; }
        }
        public IValueParameter<DoubleValue> LambdaParameter {
            get { return (IValueParameter<DoubleValue>)Parameters[LambdaParameterName]; }
        }
        public IValueParameter<IntValue> MaxNumBasisFuncsParameter {
            get { return (IValueParameter<IntValue>)Parameters[MaxNumBasisFuncsParameterName]; }
        }
        #endregion

        #region properties
        public bool ConsiderInteractions {
            get { return ConsiderInteractionsParameter.Value.Value; }
            set { ConsiderInteractionsParameter.Value.Value = value; }
        }
        public bool ConsiderDenominations {
            get { return ConsiderDenominationsParameter.Value.Value; }
            set { ConsiderDenominationsParameter.Value.Value = value; }
        }
        public bool ConsiderExponentiations {
            get { return ConsiderExponentiationsParameter.Value.Value; }
            set { ConsiderExponentiationsParameter.Value.Value = value; }
        }
        public bool ConsiderNonlinearFunctions {
            get { return ConsiderNonlinearFuncsParameter.Value.Value; }
            set { ConsiderNonlinearFuncsParameter.Value.Value = value; }
        }
        public bool ConsiderHingeFunctions {
            get { return ConsiderHingeFuncsParameter.Value.Value; }
            set { ConsiderHingeFuncsParameter.Value.Value = value; }
        }
        public double Penalty {
            get { return PenaltyParameter.Value.Value; }
            set { PenaltyParameter.Value.Value = value; }
        }
        public DoubleValue Lambda {
            get { return LambdaParameter.Value; }
            set { LambdaParameter.Value = value; }
        }
        public int MaxNumBasisFuncs {
            get { return MaxNumBasisFuncsParameter.Value.Value; }
            set { MaxNumBasisFuncsParameter.Value.Value = value; }
        }
        #endregion

        #region ctor

        [StorableConstructor]
        private FastFunctionExtraction(StorableConstructorFlag _) : base(_) { }
        public FastFunctionExtraction(FastFunctionExtraction original, Cloner cloner) : base(original, cloner) {
        }
        public FastFunctionExtraction() : base() {
            base.Problem = new RegressionProblem();
            Parameters.Add(new ValueParameter<BoolValue>(ConsiderInteractionsParameterName, "True if you want the models to include interactions, otherwise false.", new BoolValue(true)));
            Parameters.Add(new ValueParameter<BoolValue>(ConsiderDenominationParameterName, "True if you want the models to include denominations, otherwise false.", new BoolValue(true)));
            Parameters.Add(new ValueParameter<BoolValue>(ConsiderExponentiationParameterName, "True if you want the models to include exponentiation, otherwise false.", new BoolValue(true)));
            Parameters.Add(new ValueParameter<BoolValue>(ConsiderNonlinearFuncsParameterName, "True if you want the models to include nonlinear functions(abs, log,...), otherwise false.", new BoolValue(true)));
            Parameters.Add(new ValueParameter<BoolValue>(ConsiderHingeFuncsParameterName, "True if you want the models to include Hinge Functions, otherwise false.", new BoolValue(true)));
            Parameters.Add(new ValueParameter<IntValue>(MaxNumBasisFuncsParameterName, "Set how many basis functions the models can have at most. if Max Num Basis Funcs is negative => no restriction on size", new IntValue(20)));
            Parameters.Add(new OptionalValueParameter<DoubleValue>(LambdaParameterName, "Optional: the value of lambda for which to calculate an elastic-net solution. lambda == null => calculate the whole path of all lambdas"));
            Parameters.Add(new FixedValueParameter<DoubleValue>(PenaltyParameterName, "Penalty factor (alpha) for balancing between ridge (0.0) and lasso (1.0) regression", new DoubleValue(0.05)));
        }

        [StorableHook(HookType.AfterDeserialization)]
        private void AfterDeserialization() { }

        public override IDeepCloneable Clone(Cloner cloner) {
            return new FastFunctionExtraction(this, cloner);
        }
        #endregion

        public override Type ProblemType { get { return typeof(RegressionProblem); } }
        public new RegressionProblem Problem { get { return (RegressionProblem)base.Problem; } }

        protected override void Run(CancellationToken cancellationToken) {
            var models = Fit(Problem.ProblemData, Penalty, out var numBases, ConsiderExponentiations,
                ConsiderNonlinearFunctions, ConsiderInteractions,
                ConsiderDenominations, ConsiderHingeFunctions, MaxNumBasisFuncs
            );

            int i = 0;
            var numBasesArr = numBases.ToArray();
            var solutionsArr = new List<ISymbolicRegressionSolution>();

            foreach (var model in models) {
                Results.Add(new Result(
                    "Num Bases: " + numBasesArr[i++],
                    model
                ));
                solutionsArr.Add(new SymbolicRegressionSolution(model, Problem.ProblemData));
            }
            Results.Add(new Result("Model Accuracies", new ItemCollection<ISymbolicRegressionSolution>(solutionsArr)));
        }

        public static IEnumerable<ISymbolicRegressionModel> Fit(IRegressionProblemData data, double elnetPenalty, out IEnumerable<int> numBases, bool exp = true, bool nonlinFuncs = true, bool interactions = true, bool denoms = false, bool hingeFuncs = true, int maxNumBases = -1) {
            var approaches = CreateApproaches(interactions, denoms, exp, nonlinFuncs, hingeFuncs, maxNumBases, elnetPenalty);

            var allFFXModels = approaches
                .SelectMany(approach => CreateFFXModels(data, approach)).ToList();

            // Final pareto filter over all generated models from all different approaches
            var nondominatedFFXModels = NondominatedModels(data, allFFXModels);

            numBases = nondominatedFFXModels
                .Select(ffxModel => ffxModel.NumBases).ToArray();
            return nondominatedFFXModels.Select(ffxModel => ffxModel.ToSymbolicRegressionModel(data.TargetVariable));
        }

        private static IEnumerable<FFXModel> NondominatedModels(IRegressionProblemData data, IEnumerable<FFXModel> ffxModels) {
            var numBases = ffxModels.Select(ffxModel => (double)ffxModel.NumBases).ToArray();
            var errors = ffxModels.Select(ffxModel => {
                var originalValues = data.TargetVariableTestValues.ToArray();
                var estimatedValues = ffxModel.Simulate(data, data.TestIndices);
                // do not create a regressionSolution here for better performance: 
                // RegressionSolutions calculate all kinds of errors when calling the ctor, but we only need testMSE
                var testMSE = OnlineMeanSquaredErrorCalculator.Calculate(originalValues, estimatedValues, out var state);
                if (state != OnlineCalculatorError.None) throw new ArrayTypeMismatchException("could not calculate TestMSE");
                return testMSE;
            }).ToArray();

            int n = numBases.Length;
            double[][] qualities = new double[n][];
            for (int i = 0; i < n; i++) {
                qualities[i] = new double[2];
                qualities[i][0] = numBases[i];
                qualities[i][1] = errors[i];
            }

            return DominationCalculator<FFXModel>.CalculateBestParetoFront(ffxModels.ToArray(), qualities, new bool[] { false, false })
                .Select(tuple => tuple.Item1).OrderBy(ffxModel => ffxModel.NumBases);
        }

        // Build FFX models
        private static IEnumerable<FFXModel> CreateFFXModels(IRegressionProblemData data, Approach approach) {
            // FFX Step 1
            var basisFunctions = BFUtils.CreateBasisFunctions(data, approach).ToArray();

            var funcsArr = basisFunctions.ToArray();
            var elnetData = BFUtils.PrepareData(data, funcsArr);

            // FFX Step 2
            ElasticNetLinearRegression.RunElasticNetLinearRegression(elnetData, approach.ElnetPenalty, out var _, out var _, out var _, out var candidateCoeffs, out var intercept, maxVars: approach.MaxNumBases);

            // create models out of the learned coefficients
            var ffxModels = GetUniqueModelsFromCoeffs(candidateCoeffs, intercept, funcsArr, approach);
             
            // one last LS-optimization step on the training data
            foreach (var ffxModel in ffxModels) {
                if (ffxModel.NumBases > 0) ffxModel.OptimizeCoefficients(data);
            }
            return ffxModels;
        }

        // finds all models with unique combinations of basis functions
        private static IEnumerable<FFXModel> GetUniqueModelsFromCoeffs(double[,] candidateCoeffs, double[] intercept, IBasisFunction[] funcsArr, Approach approach) {
            List<FFXModel> ffxModels = new List<FFXModel>();
            List<int[]> bfCombinations = new List<int[]>();

            for (int i = 0; i < intercept.Length; i++) {
                var row = candidateCoeffs.GetRow(i);
                var nonzeroIndices = row.FindAllIndices(val => val != 0).ToArray();
                if (nonzeroIndices.Count() > approach.MaxNumBases) continue;
                // ignore duplicate models (models with same combination of basis functions)
                if (bfCombinations.Any(arr => Enumerable.SequenceEqual(arr, nonzeroIndices))) continue;
                var ffxModel = new FFXModel(intercept[i], nonzeroIndices.Select(idx => (row[idx], funcsArr[idx])));
                bfCombinations.Add(nonzeroIndices);
                ffxModels.Add(ffxModel);
            }
            return ffxModels;
        }

        private static IEnumerable<Approach> CreateApproaches(bool interactions, bool denominator, bool exponentiations, bool nonlinearFuncs, bool hingeFunctions, int maxNumBases, double penalty) {
            var approaches = new List<Approach>();
            var valids = new bool[5] { interactions, denominator, exponentiations, nonlinearFuncs, hingeFunctions };

            // return true if ALL indices of true values of arr1 also have true values in arr2
            bool follows(bool[] arr1, bool[] arr2) {
                if (arr1.Length != arr2.Length) throw new ArgumentException("invalid lengths");
                for (int i = 0; i < arr1.Length; i++) {
                    if (arr1[i] && !arr2[i]) return false;
                }
                return true;
            }

            for (int i = 0; i < 32; i++) {
                // map i to a bool array of length 5
                var arr = i.ToBoolArray(5);
                if (!follows(arr, valids)) continue;
                int sum = arr.Where(b => b).Count(); // how many features are enabled?
                if (sum >= 4) continue; // not too many features at once
                if (arr[0] && arr[2]) continue; // never need both exponent and inter
                approaches.Add(new Approach(arr[0], arr[1], arr[2], arr[3], arr[4], exponents, nonlinFuncs, maxNumBases, penalty, minHingeThr, maxHingeThr, numHingeThrs));
            }
            return approaches;
        }
    }
}