Auscitte/minimax_fit.py

## minimax_fit.py
""" Solution for the "Minimax rational fit to the exponential" problem in Stephen Boyd's Convex Optimization

    Compares a bisection-based implementation to that employing coordinate descent (iterative partial minimization).

:Copyright:
    Copyright Ry Auscitte 2021. This script is distributed under MIT License.

:Authors:
    Ry Auscitte
"""

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt


class FitProblem:
    """ An abstract class representing the rational fit to the exponential problem.

        This class implements functionality common for both methods of solving the problem.
    """
    def define_domain_bounds(self):
        self.k = 201
        self.lb = -3
        self.ub = 3

    def function_to_fit(self, x):
        return np.exp(x)

    def generate_data(self):

        self.define_domain_bounds()

        self.t = np.array([ self.lb + (self.ub - self.lb) * i / (self.k - 1.0) for i in range(self.k) ])
        self.tsq = np.multiply(self.t, self.t)

        self.y = self.function_to_fit(self.t)

    def declare_variables_and_parameters(self):

        self.a0 = cp.Variable()
        self.a1 = cp.Variable()
        self.a2 = cp.Variable()

        self.b1 = cp.Variable()
        self.b2 = cp.Variable()

        self.z = cp.Variable(self.k, pos = True)

    def declare_constraints(self):

        self.constraints = [ np.ones(self.k) + self.b1 * self.t + self.b2 * self.tsq == self.z ]

    def declare_objective(self):

        self.objective = None

        raise NotImplementedError("declare_objective() has not been implemented")

    def create_problem(self):

        self.generate_data()
        self.declare_variables_and_parameters()
        self.declare_constraints()
        self.declare_objective()
        self.prob = cp.Problem(self.objective, self.constraints)

    def call_solver(self):
        self.prob.solve(solver = cp.ECOS)

    def fitted_fun(self):
	    return np.divide(self.a0.value + self.a1.value * self.t + self.a2.value * self.tsq,
	        np.ones(self.k) + self.b1.value * self.t + self.b2.value * self.tsq)

    def evaluate_objective(self):
	    return np.linalg.norm(self.fitted_fun() - self.y, ord = np.inf)

    def print_curvature(self):
        print("Objective's curvature:", self.objective.expr.curvature)
        for i in range(len(self.constraints)):
            print("Curvature of constraint", i + 1,":", self.constraints[i].expr.curvature)


class BisectionFitProblem(FitProblem):
    """Bisection implementation
    """
    def __init__(self):
        self.eps = 0.001

    def declare_variables_and_parameters(self):

        super().declare_variables_and_parameters()
        self.s = cp.Parameter()

    def declare_objective(self):

        self.objective = cp.Minimize(0)

    def declare_constraints(self):

        super().declare_constraints()

        self.constraints.extend(
            [ self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) <= self.z * self.s,
              self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) >= -self.z * self.s,
            ])

    def solve(self, u):

        time = 0.0
        steps = 0
        l = 0.0

        while u - l > self.eps:

            self.s.value = (u + l) / 2.0

            self.call_solver()

            time += self.prob.solver_stats.solve_time
            steps += 1

            if self.prob.status == 'optimal':
                u = self.s.value
            else:
                l = self.s.value

        self.s.value = u
        self.call_solver()
        assert(self.prob.status == 'optimal')

        return (steps, time)


class BisectionFitProblemHighEps(BisectionFitProblem):
    """Bisection method with higher accuracy"""

    def __init__(self):
        self.eps = 0.0001

    def call_solver(self):
        solver = cp.ECOS
        try:
            self.prob.solve(solver = solver)
        except:
            self.prob.solve(solver = solver, abstol = 1.0e-3) #default abstol = 1.0e-7
            print("Status for the problem with numerical issues", self.prob.status)


class CoordinateDescentFitProblem(FitProblem):
    """ Impelements coordinate descent over two subsets of variables {a0, a1, a2, b1, b2} and {a0, a1, a2, s}
    """
    def __init__(self):
        self.eps = 0.0000001

    def declare_variables_and_parameters(self):

        super().declare_variables_and_parameters()

        #problem 1
        self.s_p = cp.Parameter()

        #problem 2
        self.s = cp.Variable(nonneg = True)
        self.z_p = cp.Parameter(self.k)

    def declare_objective(self):

        self.objective = cp.Minimize(0)
        self.objective2 = cp.Minimize(self.s)

    def declare_constraints(self):

        super().declare_constraints()

        self.constraints.extend(
            [ self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) <= self.s_p * self.z,
              self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) >= -self.s_p * self.z
            ])

        self.constraints2 = [ self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z_p) <= self.s * self.z_p,
                              self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z_p) >= -self.s * self.z_p
                            ]

    def create_problem(self):

        super().create_problem()
        self.prob2 = cp.Problem(self.objective2, self.constraints2)

    def call_solver(self, second = False):

        if second:
            self.prob2.solve(solver = cp.ECOS)
        else:
            self.prob.solve(solver = cp.ECOS)

    def solve(self, sv):

        self.s_p.value = sv
        val = np.finfo(float).max
        steps = 0
        time = 0.0
        objective_vals = []

        while True:

            self.call_solver()
            assert(self.prob.status == 'optimal')
            time += self.prob.solver_stats.solve_time

            self.z_p.value = np.ones(self.k) + self.b1.value * self.t + self.b2.value * self.tsq

            self.call_solver(second = True)
            assert(self.prob2.status == 'optimal')
            time += self.prob2.solver_stats.solve_time
            steps += 1

            self.s_p.value = self.prob2.value

            if abs(val - self.prob2.value) < self.eps:
                break

            val = self.prob2.value

            objective_vals.append(val)

        return (steps, time, objective_vals)


class FlawedPartialMinimizationFitProblem(FitProblem):
    """Implementation of an iterative partial minimization method (experimental).
       This method has an inherent flaw to it; it is implemented as an experiment.
    """
    def __init__(self):
        self.eps = 0.0000001
        self.steps_limit = 50

    def declare_variables_and_parameters(self):

        super().declare_variables_and_parameters()
        self.s = cp.Variable(nonneg = True)
        self.bs1 = cp.Parameter()
        self.bs2 = cp.Parameter()

    def declare_objective(self):

        self.objective = cp.Minimize(self.s)

    def declare_constraints(self):

        super().declare_constraints()

        self.constraints.extend(
            [ self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) <= self.s + self.bs1 * self.t + self.bs2 * self.tsq,
              self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) >= -(self.s + self.bs1 * self.t + self.bs2 * self.tsq)
            ])

    def compute_violations(self):
	    lhs = np.abs(self.a0.value + self.a1.value * self.t + self.a2.value * self.tsq - np.multiply(self.y, self.z.value))
	    rhs = self.s.value + self.bs1.value * self.t + self.bs2.value * self.tsq
	    return max(np.max(lhs - rhs), 0)

    def solve(self, iv1, iv2):

        self.bs1.value = iv1;
        self.bs2.value = iv2;

        val = np.finfo(float).max
        steps = 0
        time = 0.0
        objective_vals = []
        true_objective_vals = []
        violations = []

        while True:

            self.call_solver()

            time += self.prob.solver_stats.solve_time
            steps += 1

            self.bs1.value = self.s.value * self.b1.value
            self.bs2.value = self.s.value * self.b2.value

            if abs(val - self.prob.value) < self.eps:
                break

            val = self.prob.value

            objective_vals.append(val)
            violations.append(self.compute_violations())
            true_objective_vals.append(self.evaluate_objective())

            if steps > self.steps_limit:
                break

        return (steps, time, objective_vals, true_objective_vals, violations)


class FlawedPartialMinimizationFitProblemForLog(FlawedPartialMinimizationFitProblem):
    """Fits log() by experimental iterative partial minimization"""

    def define_domain_bounds(self):
        self.k = 201
        self.lb = 0.1
        self.ub = 6.1

    def function_to_fit(self, x):
        return np.log(x)


def bisection_higher_acccuracy():
    p = BisectionFitProblemHighEps()
    p.create_problem()
    steps, time = p.solve(1)
    print("Bisection method took", steps, "iterations and", time, "ms to solve the problem")
    print("Attained objective value is", p.evaluate_objective())


def main():

    #Bisection method
    print("Bisection method")

    bp = BisectionFitProblem()
    bp.create_problem()
    bp.print_curvature()
    bp_steps, bp_time = bp.solve(1) #bp.solve(0.1) runs into numerical issues; sensitive to the initial values
    print("Bisection method took", bp_steps, "iterations and", "{0:0.3f}".format(bp_time), "ms to solve the problem")
    print("Attained objective value is", bp.evaluate_objective())

    plt.plot(bp.t, bp.y, label = "target function")
    plt.plot(bp.t, bp.fitted_fun(), linestyle='dashed', label = "fitted function")
    plt.title('Bisection method')
    plt.legend()
    plt.show()

    #Coordinate descent method
    print("\nCoordinate descent method")
    cdp = CoordinateDescentFitProblem()
    cdp.create_problem()
    cdp_steps, cdp_time, cdp_obj_vals = cdp.solve(1)
    print("Coordinate descent took", cdp_steps, "iterations and", "{0:0.3f}".format(cdp_time), "ms to solve the problem")
    print("Attained objective value is", cdp.evaluate_objective())

    fig, ax = plt.subplots(1, 2)
    plt.suptitle("Coordinate descent method")
    fig.tight_layout()
    ax[0].plot(cdp_obj_vals, label = "objective (s) vs iteration number")
    ax[0].legend()
    ax[1].plot(cdp.t, cdp.y, label = "target function")
    ax[1].plot(cdp.t, cdp.fitted_fun(), linestyle = "dashed", label = "fitted function")
    ax[1].legend()
    plt.show()

    #optimal objective value is about 0.0227; setting s lower than that will produce an infeasible problem
    siv = np.concatenate((np.linspace(0.03, 0.1, 10), np.linspace(0.1, 1, 10)))
    sts = np.zeros(len(siv))
    obs = np.zeros(len(siv))
    for i in range(len(siv)):
        cdp_steps, _, _ = cdp.solve(siv[i])
        sts[i] = cdp_steps
        obs[i] = cdp.evaluate_objective()

    print("Avarage optimal objective:", "{0:0.5f}".format(np.mean(obs)), " st. deviation", "{0:0.7f}".format(np.std(obs)))

    plt.plot(siv, sts, label = "number of iterations vs initial value for s")
    plt.title('Coordinate descent method')
    plt.legend()
    plt.show()

    #Partial minimization method
    print("\nIterative partial minimization method")
    pmp = FlawedPartialMinimizationFitProblem()
    pmp.create_problem()
    pmp.print_curvature()
    pmp_steps, pmp_time, obj_vals, real_obj_vals, violations = pmp.solve(0.001, 0.001)
    print("Iterative partial minimization took", pmp_steps, "iterations and", "{0:0.3f}".format(pmp_time), "ms to solve the problem")
    print("Attained objective value is", pmp.evaluate_objective())

    fig, ax = plt.subplots(1, 2)
    plt.suptitle("Iterative partial minimization method")
    fig.tight_layout()
    ax[0].plot(obj_vals, label = "objective (s)")
    ax[0].plot(violations, label = "violations")
    ax[0].plot(real_obj_vals, label = "original objective (inf norm)")
    ax[0].legend()
    ax[1].plot(pmp.t, pmp.y, label = "target function")
    ax[1].plot(pmp.t, pmp.fitted_fun(), linestyle = "dashed", label = "fitted function")
    ax[1].legend()
    plt.show()

    #Genearing the heat map llustrating dependency between number of iterations and initial values for the problem parameters
    init_vals = np.concatenate((np.linspace(-50, -1, 10), np.linspace(-1, -0.1, 10), np.linspace(-0.1, -0.01, 10),
        np.linspace(0.01, 0.1, 10), np.linspace(0.1, 1, 10), np.linspace(1, 50, 10)))
    stephm = np.zeros((len(init_vals), len(init_vals)))
    for i in range(len(init_vals)):
        for j in range(len(init_vals)):
            st, _, _, _,_ = pmp.solve(init_vals[i], init_vals[j])
            stephm[i, j] = st

    x, y = np.meshgrid(init_vals, init_vals)

    fig, ax = plt.subplots()
    cmsh = ax.pcolormesh(x, y, stephm, cmap = 'RdYlBu_r')
    ax.set_title('Number of iterations')
    ax.axis([x.min(), x.max(), y.min(), y.max()])
    plt.ylabel('initial value of bs1')
    plt.xlabel('initial value of bs2')
    fig.colorbar(cmsh, ax = ax)
    plt.show()

    #Partial minimization method fitting log()
    print("\nIterative partial minimization method for fitting log()")
    pmplg = FlawedPartialMinimizationFitProblemForLog()
    pmplg.create_problem()
    pmplg.print_curvature()
    _, _, obj_vals, real_obj_vals, violations = pmplg.solve(0.001, 0.001)

    fig, ax = plt.subplots(1, 2)
    plt.suptitle("Iterative partial minimization method for log()")
    fig.tight_layout()
    ax[0].plot(obj_vals, label = "objective (s)")
    ax[0].plot(violations, label = "violations")
    ax[0].plot(real_obj_vals, label = "original objective (inf norm)")
    ax[0].legend()
    ax[1].plot(pmplg.t, pmplg.y, label = "target function")
    ax[1].plot(pmplg.t, pmplg.fitted_fun(), linestyle = "dashed", label = "fitted function")
    ax[1].legend()
    plt.show()


if __name__ == "__main__":
    main()
	""" Solution for the "Minimax rational fit to the exponential" problem in Stephen Boyd's Convex Optimization

	Compares a bisection-based implementation to that employing coordinate descent (iterative partial minimization).

	:Copyright:
	Copyright Ry Auscitte 2021. This script is distributed under MIT License.

	:Authors:
	Ry Auscitte
	"""

	import cvxpy as cp
	import numpy as np
	import matplotlib.pyplot as plt


	class FitProblem:
	""" An abstract class representing the rational fit to the exponential problem.

	This class implements functionality common for both methods of solving the problem.
	"""
	def define_domain_bounds(self):
	self.k = 201
	self.lb = -3
	self.ub = 3

	def function_to_fit(self, x):
	return np.exp(x)

	def generate_data(self):

	self.define_domain_bounds()

	self.t = np.array([ self.lb + (self.ub - self.lb) * i / (self.k - 1.0) for i in range(self.k) ])
	self.tsq = np.multiply(self.t, self.t)

	self.y = self.function_to_fit(self.t)

	def declare_variables_and_parameters(self):

	self.a0 = cp.Variable()
	self.a1 = cp.Variable()
	self.a2 = cp.Variable()

	self.b1 = cp.Variable()
	self.b2 = cp.Variable()

	self.z = cp.Variable(self.k, pos = True)

	def declare_constraints(self):

	self.constraints = [ np.ones(self.k) + self.b1 * self.t + self.b2 * self.tsq == self.z ]

	def declare_objective(self):

	self.objective = None

	raise NotImplementedError("declare_objective() has not been implemented")

	def create_problem(self):

	self.generate_data()
	self.declare_variables_and_parameters()
	self.declare_constraints()
	self.declare_objective()
	self.prob = cp.Problem(self.objective, self.constraints)

	def call_solver(self):
	self.prob.solve(solver = cp.ECOS)

	def fitted_fun(self):
	return np.divide(self.a0.value + self.a1.value * self.t + self.a2.value * self.tsq,
	np.ones(self.k) + self.b1.value * self.t + self.b2.value * self.tsq)

	def evaluate_objective(self):
	return np.linalg.norm(self.fitted_fun() - self.y, ord = np.inf)

	def print_curvature(self):
	print("Objective's curvature:", self.objective.expr.curvature)
	for i in range(len(self.constraints)):
	print("Curvature of constraint", i + 1,":", self.constraints[i].expr.curvature)



	class BisectionFitProblem(FitProblem):
	"""Bisection implementation
	"""
	def __init__(self):
	self.eps = 0.001

	def declare_variables_and_parameters(self):

	super().declare_variables_and_parameters()
	self.s = cp.Parameter()

	def declare_objective(self):

	self.objective = cp.Minimize(0)

	def declare_constraints(self):

	super().declare_constraints()

	self.constraints.extend(
	[ self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) <= self.z * self.s,
	self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) >= -self.z * self.s,
	])

	def solve(self, u):

	time = 0.0
	steps = 0
	l = 0.0

	while u - l > self.eps:

	self.s.value = (u + l) / 2.0

	self.call_solver()

	time += self.prob.solver_stats.solve_time
	steps += 1

	if self.prob.status == 'optimal':
	u = self.s.value
	else:
	l = self.s.value

	self.s.value = u
	self.call_solver()
	assert(self.prob.status == 'optimal')

	return (steps, time)



	class BisectionFitProblemHighEps(BisectionFitProblem):
	"""Bisection method with higher accuracy"""

	def __init__(self):
	self.eps = 0.0001

	def call_solver(self):
	solver = cp.ECOS
	try:
	self.prob.solve(solver = solver)
	except:
	self.prob.solve(solver = solver, abstol = 1.0e-3) #default abstol = 1.0e-7
	print("Status for the problem with numerical issues", self.prob.status)



	class CoordinateDescentFitProblem(FitProblem):
	""" Impelements coordinate descent over two subsets of variables {a0, a1, a2, b1, b2} and {a0, a1, a2, s}
	"""
	def __init__(self):
	self.eps = 0.0000001

	def declare_variables_and_parameters(self):

	super().declare_variables_and_parameters()

	#problem 1
	self.s_p = cp.Parameter()

	#problem 2
	self.s = cp.Variable(nonneg = True)
	self.z_p = cp.Parameter(self.k)

	def declare_objective(self):

	self.objective = cp.Minimize(0)
	self.objective2 = cp.Minimize(self.s)

	def declare_constraints(self):

	super().declare_constraints()

	self.constraints.extend(
	[ self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) <= self.s_p * self.z,
	self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) >= -self.s_p * self.z
	])

	self.constraints2 = [ self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z_p) <= self.s * self.z_p,
	self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z_p) >= -self.s * self.z_p
	]

	def create_problem(self):

	super().create_problem()
	self.prob2 = cp.Problem(self.objective2, self.constraints2)

	def call_solver(self, second = False):

	if second:
	self.prob2.solve(solver = cp.ECOS)
	else:
	self.prob.solve(solver = cp.ECOS)

	def solve(self, sv):

	self.s_p.value = sv
	val = np.finfo(float).max
	steps = 0
	time = 0.0
	objective_vals = []

	while True:

	self.call_solver()
	assert(self.prob.status == 'optimal')
	time += self.prob.solver_stats.solve_time

	self.z_p.value = np.ones(self.k) + self.b1.value * self.t + self.b2.value * self.tsq

	self.call_solver(second = True)
	assert(self.prob2.status == 'optimal')
	time += self.prob2.solver_stats.solve_time
	steps += 1

	self.s_p.value = self.prob2.value

	if abs(val - self.prob2.value) < self.eps:
	break

	val = self.prob2.value

	objective_vals.append(val)

	return (steps, time, objective_vals)



	class FlawedPartialMinimizationFitProblem(FitProblem):
	"""Implementation of an iterative partial minimization method (experimental).
	This method has an inherent flaw to it; it is implemented as an experiment.
	"""
	def __init__(self):
	self.eps = 0.0000001
	self.steps_limit = 50

	def declare_variables_and_parameters(self):

	super().declare_variables_and_parameters()
	self.s = cp.Variable(nonneg = True)
	self.bs1 = cp.Parameter()
	self.bs2 = cp.Parameter()

	def declare_objective(self):

	self.objective = cp.Minimize(self.s)

	def declare_constraints(self):

	super().declare_constraints()

	self.constraints.extend(
	[ self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) <= self.s + self.bs1 * self.t + self.bs2 * self.tsq,
	self.a0 + self.a1 * self.t + self.a2 * self.tsq - cp.multiply(self.y, self.z) >= -(self.s + self.bs1 * self.t + self.bs2 * self.tsq)
	])

	def compute_violations(self):
	lhs = np.abs(self.a0.value + self.a1.value * self.t + self.a2.value * self.tsq - np.multiply(self.y, self.z.value))
	rhs = self.s.value + self.bs1.value * self.t + self.bs2.value * self.tsq
	return max(np.max(lhs - rhs), 0)

	def solve(self, iv1, iv2):

	self.bs1.value = iv1;
	self.bs2.value = iv2;

	val = np.finfo(float).max
	steps = 0
	time = 0.0
	objective_vals = []
	true_objective_vals = []
	violations = []

	while True:

	self.call_solver()

	time += self.prob.solver_stats.solve_time
	steps += 1

	self.bs1.value = self.s.value * self.b1.value
	self.bs2.value = self.s.value * self.b2.value

	if abs(val - self.prob.value) < self.eps:
	break

	val = self.prob.value

	objective_vals.append(val)
	violations.append(self.compute_violations())
	true_objective_vals.append(self.evaluate_objective())

	if steps > self.steps_limit:
	break

	return (steps, time, objective_vals, true_objective_vals, violations)



	class FlawedPartialMinimizationFitProblemForLog(FlawedPartialMinimizationFitProblem):
	"""Fits log() by experimental iterative partial minimization"""

	def define_domain_bounds(self):
	self.k = 201
	self.lb = 0.1
	self.ub = 6.1

	def function_to_fit(self, x):
	return np.log(x)



	def bisection_higher_acccuracy():
	p = BisectionFitProblemHighEps()
	p.create_problem()
	steps, time = p.solve(1)
	print("Bisection method took", steps, "iterations and", time, "ms to solve the problem")
	print("Attained objective value is", p.evaluate_objective())



	def main():

	#Bisection method
	print("Bisection method")

	bp = BisectionFitProblem()
	bp.create_problem()
	bp.print_curvature()
	bp_steps, bp_time = bp.solve(1) #bp.solve(0.1) runs into numerical issues; sensitive to the initial values
	print("Bisection method took", bp_steps, "iterations and", "{0:0.3f}".format(bp_time), "ms to solve the problem")
	print("Attained objective value is", bp.evaluate_objective())

	plt.plot(bp.t, bp.y, label = "target function")
	plt.plot(bp.t, bp.fitted_fun(), linestyle='dashed', label = "fitted function")
	plt.title('Bisection method')
	plt.legend()
	plt.show()

	#Coordinate descent method
	print("\nCoordinate descent method")
	cdp = CoordinateDescentFitProblem()
	cdp.create_problem()
	cdp_steps, cdp_time, cdp_obj_vals = cdp.solve(1)
	print("Coordinate descent took", cdp_steps, "iterations and", "{0:0.3f}".format(cdp_time), "ms to solve the problem")
	print("Attained objective value is", cdp.evaluate_objective())

	fig, ax = plt.subplots(1, 2)
	plt.suptitle("Coordinate descent method")
	fig.tight_layout()
	ax[0].plot(cdp_obj_vals, label = "objective (s) vs iteration number")
	ax[0].legend()
	ax[1].plot(cdp.t, cdp.y, label = "target function")
	ax[1].plot(cdp.t, cdp.fitted_fun(), linestyle = "dashed", label = "fitted function")
	ax[1].legend()
	plt.show()

	#optimal objective value is about 0.0227; setting s lower than that will produce an infeasible problem
	siv = np.concatenate((np.linspace(0.03, 0.1, 10), np.linspace(0.1, 1, 10)))
	sts = np.zeros(len(siv))
	obs = np.zeros(len(siv))
	for i in range(len(siv)):
	cdp_steps, _, _ = cdp.solve(siv[i])
	sts[i] = cdp_steps
	obs[i] = cdp.evaluate_objective()

	print("Avarage optimal objective:", "{0:0.5f}".format(np.mean(obs)), " st. deviation", "{0:0.7f}".format(np.std(obs)))

	plt.plot(siv, sts, label = "number of iterations vs initial value for s")
	plt.title('Coordinate descent method')
	plt.legend()
	plt.show()

	#Partial minimization method
	print("\nIterative partial minimization method")
	pmp = FlawedPartialMinimizationFitProblem()
	pmp.create_problem()
	pmp.print_curvature()
	pmp_steps, pmp_time, obj_vals, real_obj_vals, violations = pmp.solve(0.001, 0.001)
	print("Iterative partial minimization took", pmp_steps, "iterations and", "{0:0.3f}".format(pmp_time), "ms to solve the problem")
	print("Attained objective value is", pmp.evaluate_objective())

	fig, ax = plt.subplots(1, 2)
	plt.suptitle("Iterative partial minimization method")
	fig.tight_layout()
	ax[0].plot(obj_vals, label = "objective (s)")
	ax[0].plot(violations, label = "violations")
	ax[0].plot(real_obj_vals, label = "original objective (inf norm)")
	ax[0].legend()
	ax[1].plot(pmp.t, pmp.y, label = "target function")
	ax[1].plot(pmp.t, pmp.fitted_fun(), linestyle = "dashed", label = "fitted function")
	ax[1].legend()
	plt.show()

	#Genearing the heat map llustrating dependency between number of iterations and initial values for the problem parameters
	init_vals = np.concatenate((np.linspace(-50, -1, 10), np.linspace(-1, -0.1, 10), np.linspace(-0.1, -0.01, 10),
	np.linspace(0.01, 0.1, 10), np.linspace(0.1, 1, 10), np.linspace(1, 50, 10)))
	stephm = np.zeros((len(init_vals), len(init_vals)))
	for i in range(len(init_vals)):
	for j in range(len(init_vals)):
	st, _, _, _,_ = pmp.solve(init_vals[i], init_vals[j])
	stephm[i, j] = st

	x, y = np.meshgrid(init_vals, init_vals)

	fig, ax = plt.subplots()
	cmsh = ax.pcolormesh(x, y, stephm, cmap = 'RdYlBu_r')
	ax.set_title('Number of iterations')
	ax.axis([x.min(), x.max(), y.min(), y.max()])
	plt.ylabel('initial value of bs1')
	plt.xlabel('initial value of bs2')
	fig.colorbar(cmsh, ax = ax)
	plt.show()

	#Partial minimization method fitting log()
	print("\nIterative partial minimization method for fitting log()")
	pmplg = FlawedPartialMinimizationFitProblemForLog()
	pmplg.create_problem()
	pmplg.print_curvature()
	_, _, obj_vals, real_obj_vals, violations = pmplg.solve(0.001, 0.001)

	fig, ax = plt.subplots(1, 2)
	plt.suptitle("Iterative partial minimization method for log()")
	fig.tight_layout()
	ax[0].plot(obj_vals, label = "objective (s)")
	ax[0].plot(violations, label = "violations")
	ax[0].plot(real_obj_vals, label = "original objective (inf norm)")
	ax[0].legend()
	ax[1].plot(pmplg.t, pmplg.y, label = "target function")
	ax[1].plot(pmplg.t, pmplg.fitted_fun(), linestyle = "dashed", label = "fitted function")
	ax[1].legend()
	plt.show()


	if __name__ == "__main__":
	main()