/*demo_flutter_density.cpp*/
/*
* Copyright 2014 IPOL Image Processing On Line http://www.ipol.im/
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

/**
 * @mainpage Flutter shutter density estimator
 *
 * README.txt:
 * @verbinclude README.txt
 */

/**
 * @file demo_flutter_density.cpp
 * @brief Flutter shutter density estimator
 * @param argv[1]: char containg the file name of the text file containing
 * the flutter ode coefficients \alpha_k
 * @param argv[2]:  char containing the name of the text file to store
 * the estimated probability density associated to the code.
 * @param argv[3]: flag log scale 1 to get log(1+\rho(v))
 * @author Yohann Tendero <tendero@math.ucla.edu>
 */

#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>
#include "routines.h"
using namespace std;
#ifndef M_PI
/**
 * M_PI is a POSIX definition for Pi
 */
#define M_PI 3.14159265358979323846
#endif


/**
 * @brief main function call
 *
INPUTs:
\Delta t, and the file name containing the flutter code, flag for Log(1+\rho(v))

output_rho : file text.
 *
 */

int main(int argc, char **argv)
{
//Check arguments : IN OUT;
    if (argc < 3)
    {
    std::cerr << " ******************************************* " << std::endl
    << " **********  Flutter shutter density estimator  ****** " << std::endl
    << " EX: ./flutter_estimator snapshot.txt output_rho.txt 0 " << std::endl
    << " ***************************************************** " << std::endl
    << " *************  Yohann Tendero, 2014  **************** " << std::endl
    << " ***************************************************** " << std::endl;
        return 1;
    }


    /**
    * @brief : Overall description:
    * Read code, compute |\hat \alpha|, differentiate. Normalize.
    *
    **/




    ////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////
    //// Step1. Reading arguments, init ////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////

    char* flutter_code_file=argv[1]; // file_name containing the flutter code
    char* output_rho_file=argv[2];
    vector<double> code; // to store the flutter code
    double num_plot_points=1001; // number of points for the plot
    vector<double> output_rho;  // to store the estimated density
    vector<double> output_wprime;  // to store the estimated w'
    double epsilon=2/num_plot_points; // precision for numerical
    //differentiation estimation

    ////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////
    //// Step2 Read file put flutter code in vector  ///////////////////////////
    ////////////////////////////////////////////////////////////////////////////

    code=fileToVector(flutter_code_file);

    ////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////
    //// Step3 density estimation  /////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////
    double probability_of_velocity; // to store $\rho(velocity)$

    /// The following loop implements Algo 2, step 4.
    for(double velocity=-1; velocity<=1;
            velocity=velocity+2/num_plot_points)
    {
        probability_of_velocity=-velocity*w_prime_estimator(M_PI*velocity,
                                code, epsilon); // this is Thm 2 eq (26)
        if(probability_of_velocity>=0)//remove negative values
            output_rho.push_back(probability_of_velocity);
        else output_rho.push_back(0); // this implements the variant given
                                       // by eq (27).
    }
    // renormalization so \int \rho(v) dv =1
    double sum=0;
    for(unsigned k=0; k<output_rho.size(); k++)
        sum=sum+output_rho[k];
    sum=sum*2/num_plot_points; // Length of the interval/ Number of points.

    for(unsigned k=0; k<output_rho.size(); k++)
        output_rho[k]=output_rho[k]/sum;
    if(atof(argv[3])==1)
    {
        for(unsigned k=0; k<output_rho.size(); k++)
            output_rho[k]=log(output_rho[k]+1); //for log scaling the output_rho
    }
    ////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////
    //// Step4 : Writing the outputs  //////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////
    double w_prime_velocity; // to store $w'(velocity)$

    for(double velocity=0; velocity<=1; velocity=velocity+2/num_plot_points)
    {
        w_prime_velocity=w_prime_estimator(M_PI*velocity, code, epsilon);
        output_wprime.push_back(w_prime_velocity);

    }

    ///WRITE the output_rho (the density) to file : output_rho_file
    ofstream file_output_rho(output_rho_file, ios::out | ios::trunc);
    // opening the file (erase if previous exists)
    if (file_output_rho)
    {
        int k=0;
        for(double velocity=-1; velocity<=1;
                velocity=velocity+2/num_plot_points)
        {
            file_output_rho << velocity  <<" " << output_rho[k] << endl;
            k++;
        }

        file_output_rho.close();
    }
    else
        cerr << "Unable to open file !" << endl;


    ///WRITE the estimated w' to file : "wprime.txt"
    ofstream file_output_wprime("wprime.txt", ios::out | ios::trunc);
    // opening the file (erase if previous exists)
    if ("wprime.txt")
    {
        int k=0;
        for(double velocity=0; velocity<=1; velocity=velocity+2/num_plot_points)
        {
            file_output_wprime << velocity  <<" " << output_wprime[k] << endl;
            k++;
        }

        file_output_wprime.close();
    }
    else
        cerr << "Unable to open file !" << endl;



    return 0;
}