/*****************************************************************************
 *                                                                          
 * This software module was originally developed by 
 *
 * J.I. Ronda (UPM-GTI / ACTS-MoMuSys)
 *
 * Copyright (c) 1997
 *
 *****************************************************************************/

/***********************************************************HeaderBegin*******
 *                                                                         
 * File:        rc_upm_model2.c
 *
 * Author:      J.I. Ronda, UPM-GTI
 * Created:     05.12.97
 *                                                                         
 * Description: UPM prediction model functions
 *
 * Notes:  These functions make use of NUMERICAL RECIPIES     
 *
 * Flags:       -D_RC_DEBUG_  -  RC debugging   
 *
 * Modified:
 *
 ***********************************************************HeaderEnd*********/

/************************    INCLUDE FILES    ********************************/
/*#ifdef _RC_UPM2_*/

#include "rc.h"
#include "momusys.h"
#include "mom_vop.h"
#include "vm_config.h"
#include "mom_vol.h"
#include "mom_structs.h"
#include "vm_enc_defs.h"
#include "rc_util2.h"


/***********************************************************CommentBegin******
 * 
 * -- rc_comp_cte --
 *
 * Author : 
 *      J. Ignacio Ronda, UPM-GTI     
 *
 * Created :            
 *      09-12-97
 *
 * Purpose :
 *      Computes the optimum cte for model
 *           y[i] = p1*x1[i] + p2*x2[i] + cte
 *
 ***********************************************************CommentEnd********/
double rc_comp_cte(double *y, double *x1, double *x2,
                   double p1, double p2, int n)
{
double sum=0, p3, perr, average;
int i;

for (i=0; i<n; i++)
   sum += (y[i]-p1*x1[i]-p2*x2[i]);

p3 = sum/n;

perr = 0;
average = 0;
for (i=0;i<n;i++)
  {
  Double pred;
  pred = p1*x1[i]+p2*x2[i]+p3;
  perr += CUAD(y[i]-pred);
  average += y[i];
  }
perr /= n;
average /= n;
/*JIR*/fprintf(stdout, "JIR: rel. perr (2) = %E\n", sqrt(perr)/average);

return p3;

}

/***********************************************************CommentBegin******
 *
 * -- rc_update_UPM_model2 --
 *
 * Author : 
 *      Martina Eckert, UPM-GTI     
 *
 * Created :            
 *      18-06-97
 *
 * Purpose :
 *      Update of UPM model. Different for rate and distortion. 
 * 
 * Arguments in :       
 *      RC_HIST *rch    Static history struct
 *      Int     type    rate or distortion parameters 
 *
 * Arguments in/out :   
 *      Double  *p1 <-> Params for prediction
 *      Double  *p2      
 *      Double  *p3      
 *
 * Arguments out :      
 *      -
 *
 * Return values :      
 *      -
 *
 * Side effects :       
 *      -
 *
 * Description :        
 *
 *
 * See also :
 *
 *
 * Modified : 3.12.97 M. Eckert: Include model improvement from JIR:
 *                    Calculation of p1,p2,p3 with SVD, if monotony criterion
 *                    is not fulfulled, adaption to previous parameters
 *
 *
 ***********************************************************CommentEnd********/

 void rc_update_UPM_model2(
			   RC_HIST     *rch,
			   Double      *p1,   
			   Double      *p2,   
			   Double      *p3,   
			   Int         type   
			   )
{

   Double E[MAX_SLIDING_WINDOW+10];   /* R/E or D/E data */
   Double QP[MAX_SLIDING_WINDOW+10];  /* Rate pred.: inverse QP, dist. pred.: QP */
   Double QP2[MAX_SLIDING_WINDOW+10]; /* Rate pred.: inverse QP^2, dist. pred.: QP^2 */
   Double ONE[MAX_SLIDING_WINDOW+10]; /* Rate pred.: Ones */
   Double *data[3];

   Int    n;                   /* Number of data to employ */
   Int    i, count,  model_ok;
   Double x1=99., x2=99./*, x3=99.*/;
   Double x[3];
   /*   Int range;*/
   Double alpha, alpha_inc=1./100.;
   Double pp1, pp2;
   Double x_min = 0, x_max = 0;
   Int    *bits_text;
   Double *dist;
   Double *mad_text;
   Int    *num_pel;
   Int    *qp;


#ifdef _RC_DEBUG_
   Int     last_qp;
   /*** JIR ***/static Int mcount=0;

   fprintf(stdout, "RC: RC_UPM: --> rc_update_UPM_model2\n");
   fflush(stdout);
   rch_print(rch, stdout, "rc_update_UPM_model2");

   /*** JIR ***/if (type == 1) mcount++;
   /***/printf("JIR%d<%d>Y: p1= %f ; p2= %f ; p3= %f;\n", type, mcount, *p1, *p2, *p3);
#endif

   /* Computation of number of data to employ */

   n = MIN(rch_get_n(rch), MAX_SLIDING_WINDOW);

   /* Array initialization */

   bits_text = rch_get_bits_text(rch);
   dist      = rch_get_dist(rch);
   mad_text  = rch_get_mad_text(rch);
   num_pel   = rch_get_pixels(rch);
   qp        = rch_get_qp(rch);

#ifdef _RC_DEBUG_
   last_qp   = rch_get_last_qp(rch);
  /***/printf("JIR%d<%d>: data={\n", type, mcount);
#endif

   for (i = rch_dec_mod(rch,rch_get_ptr(rch)), count = 0; count < n; count++, 
           i = rch_dec_mod(rch,i))
   {
      switch (type) 
         {
         case 1:
            E[count]  = (Double)bits_text[i] / ( mad_text[i] * num_pel[i] ); 
            QP[count]  = 1./(Double)qp[i];
            QP2[count] = QP[count] * QP[count];
            ONE[count] = 1.;
            break;
         case 2:
            E[count]  = 255./pow(10., dist[i] / 20.);
            QP[count]  = (Double)qp[i];
            QP2[count] = QP[count] * QP[count];
            ONE[count] = 1.;
            break;
         default:
            error_exit("RC_update_UPM_model2: Error: Mode not supported\n");
         }
#ifdef _RC_DEBUG_
         /***/if (count < n-1)
         /***/   printf("JIR%d<%d>: {%f,%f},\n", type, mcount, QP[count], E[count]);
         /***/else
         /***/   printf("JIR%d<%d>: {%f,%f}}\n", type, mcount, QP[count], E[count]);
#endif
      
   }
   /* JIR: Repetition of last data */
   E[count] = E[0];
   QP[count] = QP[0];
   QP2[count] = QP2[0];
   ONE[count] = ONE[0];
   count++;
   E[count] = E[0];
   QP[count] = QP[0];
   QP2[count] = QP2[0];
   ONE[count] = ONE[0];

#ifdef _RC_DEBUG_
   fprintf(stdout, "Computation of prediction parameters (Sarnoff) \n");
   
   /* LS computation */

   /***/printf("--> rc_nsvd_comp\n");
   /***/fflush(stdout);
#endif

   switch (n)
      {
      case 1:
         switch (type)
            {
            case 1: /* Rate, one data */
               data[0] = QP;
               /***
               rc_nsvd_comp(E, data, n, 1, x);
               x1 = x[0];
               ***/
               x1 = E[0]/QP[0];
               x2 = 0;
               /*x3 = 0;*/
               break;
           case 2: /* Dist., one data */
               data[0] = QP2;
               /***
               rc_nsvd_comp(E, data, n, 1, x);
               x2 = x[0];
               ***/
               x1 = 0;
               x2 = E[0]/QP2[0];
               /*x3 = 0;*/
               break;
           default:
               error_exit("RC_update_UPM_model2: Error: Mode not supported\n");
           }
         break;
      default:
         data[0] = QP;
         data[1] = QP2;
         data[2] = ONE;
         rc_nsvd_comp(E, data, n+2, 3, x);
         x1 = x[0];
         x2 = x[1];
         /*x3 = x[2];*/
      }

#ifdef _RC_DEBUG_
   /***/printf("<-- rc_nsvd_comp\n");
   /***/fflush(stdout);
   /*** printf("JIR%d<%d>X: x1= %f ; x2= %f ; x3= %f;\n", type, mcount, x1, x2, x3);*/

   fprintf(stdout, "RC: RC_Q2: <-- rc_update_2ls_comp: p1= %f  p2= %f \n", 
           x1, x2);
#endif

   /* Outlier removing (arrays reinitialization) */

   /*** JIR ***
   rc_2ls_remove_out(E, QP, QP2, &n, &x1, &x2);

   * LS recomputation *

   rc_2ls_comp_upm(E, QP, QP2, n, &x1, &x2, &range);
   *** JIR ***/

#ifdef _RC_DEBUG_
   fprintf(stdout, "RC: RC_Q2: <-- rc_update_2ls_comp: p1= %f  p2= %f \n", 
           x1, x2);
#endif
       
   switch (type)
      {
      case 1:
        /***
        x_min = 1./MAX(QP[0]-QP_INC, MIN_QUANT);
        x_max = 1./MIN(QP[0]+QP_INC, MAX_QUANT);
        ***/
        x_min = 1./MIN_QUANT;
        x_max = 1./MAX_QUANT;
        break;
      case 2:
        /***
        x_min = MAX(QP[0]-QP_INC, MIN_QUANT);
        x_max = MIN(QP[0]+QP_INC, MAX_QUANT);
        ***/
        x_min = MIN_QUANT;
        x_max = MAX_QUANT;
        break;
      default:
        error_exit("RC_update_UPM_model2: Error: Mode not supported\n");
      }

   model_ok = 0;
   for (alpha=0; alpha<=1 && !model_ok; alpha+=alpha_inc)
      {
      pp1 = alpha*(*p1)+(1-alpha)*x1;
      pp2 = alpha*(*p2)+(1-alpha)*x2;
      model_ok = rc_decide(pp1, pp2, x_min, x_max);
      }
#ifdef _RC_DEBUG_
   /*** JIR ***/printf("JIR%d<%d>: alpha= %f\n", type, mcount, alpha);
#endif
   *p3 = rc_comp_cte(E, QP, QP2, pp1, pp2, n);
 
   /*** JIR
   if (!model_ok)
   {
      if (n > 1)
      {
         rc_linear(E, QP, QP2, n, &x1, &x2, &x3);
#ifdef _RC_DEBUG_
         fprintf(stdout, "RC: RC_Q2: <-- rc_update_2ls_comp: p1= %f, p2= %f, p3= %f \n", 
                 x1, x2, x3);
#endif
      }
      if (n <= 1 || x1 <= 0)
      {
#ifdef _RC_DEBUG_
         fprintf(stdout,"RC: Model not changed! \n"); fflush(stdout);
         fprintf(stdout, "RC: RC_Q2: <-- rc_update_2ls_comp: p1= %f, p2= %f, p3= %f \n", 
	 *p1, *p2, *p3);
#endif
         return;
      }
   }
   ***/

   *p1 = pp1;
   *p2 = pp2;

#ifdef _RC_DEBUG_

  /***/printf("JIR%d<%d>Y: qp= %d  y1= %f ; y2= %f ; y3= %f;\n", type, mcount, last_qp, *p1, *p2, *p3);
   
   fprintf(stdout, "RC: RC_Q2: <-- rc_update_2ls_comp: p1= %f, p2= %f, p3= %f \n", 
           *p1, *p2, *p3);
#endif

}
/*#endif*/

