/*
 * Copyright 2011-2013 JOptimizer
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */
package com.joptimizer.optimizers;

import org.apache.commons.lang3.ArrayUtils;
import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;

import cern.colt.matrix.DoubleFactory1D;
import cern.colt.matrix.DoubleFactory2D;
import cern.colt.matrix.DoubleMatrix1D;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.linalg.Algebra;

import com.joptimizer.functions.ConvexMultivariateRealFunction;
import com.joptimizer.functions.FunctionsUtils;
import com.joptimizer.functions.LinearMultivariateRealFunction;
import com.joptimizer.util.Utils;

/**
 * Basic Phase I Method (implemented as a Primal-Dual Method).
 * 
 * @see "S.Boyd and L.Vandenberghe, Convex Optimization, p. 579"
 * @author alberto trivellato (alberto.trivellato@gmail.com)
 */
public class BasicPhaseIPDM {
	
	private PrimalDualMethod originalProblem;
	private int originalDim =-1;
	private int dim =-1;
	private Algebra ALG = Algebra.DEFAULT;
	private DoubleFactory1D F1 = DoubleFactory1D.dense;
	private DoubleFactory2D F2 = DoubleFactory2D.dense;
	private Log log = LogFactory.getLog(this.getClass().getName());

	public BasicPhaseIPDM(PrimalDualMethod originalProblem) {
		this.originalProblem = originalProblem;
		originalDim = originalProblem.getDim();
		this.dim = originalProblem.getDim()+1;//variable Y=(X, s)
	}
	
	public DoubleMatrix1D findFeasibleInitialPoint() throws Exception{
		log.debug("findFeasibleInitialPoint");
		
		OptimizationRequest or = new OptimizationRequest();
		
		//objective function: s
		DoubleMatrix1D C = F1.make(dim);
		C.set(dim-1, 1.);
		LinearMultivariateRealFunction objectiveFunction = new LinearMultivariateRealFunction(C.toArray(), 0);
		or.setF0(objectiveFunction);
		or.setToleranceFeas(originalProblem.getToleranceFeas());
		or.setTolerance(originalProblem.getTolerance());
		
	  // Inquality constraints: fi(X)-s
		ConvexMultivariateRealFunction[] inequalities = new ConvexMultivariateRealFunction[originalProblem.getFi().length];
		for(int i=0; i originalProblem.getToleranceFeas()){
				throw new Exception("The initial point for Basic Phase I Method must be equalities-feasible");
			}
		}
		
		DoubleMatrix1D originalFiX0 = originalProblem.getFi(X0);
		
		//lucky strike?
		int maxIneqIndex = Utils.getMaxIndex(originalFiX0.toArray());
		if(originalFiX0.get(maxIneqIndex) + originalProblem.getTolerance()<0){
			//the given notFeasible starting point is in fact already feasible
			return X0;
		}
		
		//DoubleMatrix1D initialPoint = F1.make(1, -Double.MAX_VALUE);
		DoubleMatrix1D initialPoint = F1.make(1, Math.sqrt(originalProblem.getToleranceFeas()));
		initialPoint = F1.append(X0, initialPoint);
		for(int i=0; i0){
		if(ineq.get(maxIneqIndex)>=0){	
			throw new Exception("Infeasible problem");
		}

        return ret;
	}
	
	private class PhaseIPrimalDualMethod extends PrimalDualMethod{
		@Override
		protected boolean checkCustomExitConditions(DoubleMatrix1D Y){
			DoubleMatrix1D X = Y.viewPart(0, dim-1);
			DoubleMatrix1D ineqX = originalProblem.getFi(X);
			int ineqMaxIndex = Utils.getMaxIndex(ineqX.toArray());
			boolean b1 = (ineqX.get(ineqMaxIndex) + getTolerance() <0) || Y.get(Y.size()-1)<0;
			DoubleMatrix1D originalRPriX = F1.make(0);
			if(getA()!=null){
				originalRPriX = originalProblem.getA().zMult(X, originalProblem.getB().copy(), 1., -1., false);
			}
			boolean b2 = Math.sqrt(ALG.norm2(originalRPriX)) < originalProblem.getToleranceFeas();
			log.debug("b1: " + b1);
			log.debug("b2: " + b2);
			log.debug("checkCustomExitConditions: " + (b1 && b2));
			return b1 && b2;
		}
	}
	
	/**
	 * Just looking for one out of all the possible solutions.
	 * @see "Convex Optimization, C.5 p. 681".
	 */
	private DoubleMatrix1D findOneRoot(double[][] A, double[] b) throws Exception{
		return F1.make(originalProblem.findEqFeasiblePoint(A, b));
	}
}