/**
 * - méthodes de pénalités intérieures
 * - méthodes de pénalités extérieures
 * fonction objectif : R -> R : f(x) = 2x
 *      * contrainte x >= 1.
 *
 * Pour la méthode des pénalités extérieures :
 * - fonction de pénalité (max(0,1-x))^2.
 */

#include <stdio.h>
#include <math.h>

#include "penality.h"

#define MAX(a,b) (a<b)?b:a;

typedef struct _vector
{
    int      dimension;
    double * components;
} Vector;

static double solution = 1;

double error(double value)
{
    return fabs(solution - value);
}

double power(double x, int pow)
{
    int i;
    double result = x;
    for(i=2;i<=pow;i++)
        result *= x;
    return result;
}

double j(double x)
{
    return 2 * x;
}

double g(double x, double epsilon)
{
    return j(x) + ((1.0/epsilon) * (1.0/(1.0 - x)));
}

double g_1(double x, double epsilon)
{
    return 2.0 + (1.0/epsilon) * (x/power(1-x,2));
}

double resolveGradientPasFixe(double x0[], double pas0, double espilon, int max_iter, double epsilon_penality)
{
    double current_iteration, old_iteration;
    double norme, diff;
    float old_error = -1, current_error;
    int counter = 0,i;

    printf("\nVecteur initial x : %lf\n",x0[0]);

    old_iteration = x0[0];

    printf("Init old iteration : %f\n",old_iteration);

    do
    {
        current_iteration = old_iteration - (pas0 * g_1(old_iteration,epsilon_penality));


        diff = fabs(current_iteration - old_iteration);

        //Calcul de l'erreur
        current_error = error(current_iteration);

        // Info debug
        printf("\n [Current iteration] : %lf -> %lf",current_iteration,g(current_iteration,epsilon_penality));
        printf("\n [Old iteration] : %lf -> %lf",old_iteration,g(old_iteration,epsilon_penality));
        printf("\n \tErreur : %f",current_error);
        if(old_error != -1)
            printf(" [ difference ( %f ) ]",diff);

        // Préparation prochaine itération
        old_iteration = current_iteration;

        // Actualisation statistique sur l'erreur
        old_error = current_error;

    }
    while( diff > espilon && ++counter < max_iter );

    // Info Debug
    printf("\nCounter value : %d",counter);
    printf("\nDifference : %f",diff);
    printf("\nVecteur optimum : %lf",current_iteration);
    printf("\nFunction value for x : %f",g(current_iteration,epsilon_penality));

    return current_iteration;
}

double resolveByPenality(int dimension, double x[], double epsilon, int max_iter)
{
    Vector iteration; // vecteur
    double current_result, old_result;       // resultat

    iteration.dimension = dimension;
    iteration.components = x;

    double epsilon_penality = 20, diff = 100;
    int counter = 1;

    do
    {
        printf("Iteration %d\n",counter);
        printf("Epsilon actuel : %lf\n",epsilon_penality);
        printf("Valeur actuelle : %lf\n",iteration.components);

        current_result = resolveGradientPasFixe(iteration.components,0.3,EPSILON,MAX_ITER,epsilon_penality);

        printf("\tResultat du sous-problème : %lf\n",current_result);

        if(old_result != -1)
            diff = fabs(current_result - old_result);

        printf("\tDifférence : %lf\n",diff);
        
        old_result = current_result;

        epsilon_penality /= 2;

    }while(counter++ < max_iter && diff > epsilon);


    printf("\nFin\n");
    printf("Nombre d'itération : %ld\n",counter);
    printf("Optimum : %ld\n",current_result);
    printf("Difference : %lf\n",diff);
    printf("Epsilon : %lf\n",epsilon_penality);

    return current_result;
}