df1b2lp9.cpp

Go to the documentation of this file.

/*
 * $Id$
 *
 * Author: David Fournier
 * Copyright (c) 2008-2012 Regents of the University of California
 */
#include <admodel.h>
#include <df1b2fun.h>
#include <adrndeff.h>
#ifdef DEBUG
  #include <cassert>
#endif

static int no_stuff=0;
/*
static void print(double ff,dvector& uuu,dvector& gg)
{
  cout << setprecision(10) << setw(19) << ff << " "
       << setw(19) << uuu   << "  "  << setw(19) << gg << endl;
}
*/

typedef fmm* pfmm;

dvector laplace_approximation_calculator::get_uhat_quasi_newton_block_diagonal
  (const dvector& x,function_minimizer * pfmin)
{
  if (separable_function_difference)
  {
    delete separable_function_difference;
    separable_function_difference=0;
  }
#ifdef DEBUG
  assert(num_separable_calls > 0);
#endif

  separable_function_difference = new dvector(1,num_separable_calls);

  fmm** pfmc1 = new pfmm[static_cast<unsigned int>(num_separable_calls)];
  pfmc1--;
  ivector ishape(1,num_separable_calls);
  dvector gmax(1,num_separable_calls);
  gmax.initialize();

  int* pishapei = ishape.get_v() + 1;
  for (int i=1;i<=num_separable_calls;i++)
  {
    int m = (*derindex)(i).indexmax();
    *pishapei = m;
    if (m>0)
    {
      fmm* pfmc1i = new fmm(m);
      pfmc1i->iprint=0;
      pfmc1i->crit=inner_crit;
      pfmc1i->ireturn=0;
      pfmc1i->itn=0;
      pfmc1i->ifn=0;
      pfmc1i->ialph=0;
      pfmc1i->ihang=0;
      pfmc1i->ihflag=0;
      pfmc1i->maxfn=100;
      pfmc1i->gmax=1.e+100;
      pfmc1i->use_control_c=0;
      pfmc1[i] = pfmc1i;
    }
    else
    {
      pfmc1[i]= (fmm *)(0);
    }

    ++pishapei;
  }
  dmatrix gg(1,num_separable_calls,1,ishape);
  dmatrix ggb(1,num_separable_calls,1,ishape);
  dmatrix uu(1,num_separable_calls,1,ishape);
  dmatrix uub(1,num_separable_calls,1,ishape);
  dvector ff(1,num_separable_calls);
  dvector ffb(1,num_separable_calls);
  ivector icon(1,num_separable_calls);
  icon.initialize();
  ffb=1.e+100;

  double f=0.0;
  double fb=1.e+100;
  dvector g(1,usize);
  dvector ub(1,usize);
  independent_variables u(1,usize);
  gradcalc(0,g);
  fmc1.itn=0;
  fmc1.ifn=0;
  fmc1.ireturn=0;
  initial_params::xinit(u);    // get the initial values into the
  fmc1.ialph=0;
  fmc1.ihang=0;
  fmc1.ihflag=0;

  if (init_switch)
  {
    u.initialize();
  }

  std::ostream& output_stream = get_output_stream();
  for (int ii=1;ii<=2;ii++)
  {
    // get the initial u into the uu's
    dvector* puui = &uu(1);
    for (int i=1;i<=num_separable_calls;i++)
    {
      int m=(*derindex)(i).indexmax();
      double* puuij = puui->get_v() + 1;
      for (int j=1;j<=m;j++)
      {
        *puuij = u((*derindex)(i)(j));

        ++puuij;
      }
      ++puui;
    }

#ifdef DIAG
    bool loop_flag = false;
    int loop_counter = 0;
#endif

    fmc1.dfn=1.e-2;
    dvariable pen=0.0;
    int converged=0;
    int initrun_flag=1;
    while (converged==0)
    {
#ifdef DIAG
      if (loop_flag) loop_counter++;
      if (loop_counter > 18)
      {
        cout << loop_counter;
      }
#endif
      if (!initrun_flag)
      {
        converged=1;
      }
      for (int i=1;i<=num_separable_calls;i++)
      {
        if (ishape(i)>0) //check to see if there are any active randoem effects
        {               // in this function call
          if (!icon(i))
          {
            independent_variables& uuu=*(independent_variables*)(&(uu(i)));
            (pfmc1[i])->fmin(ff[i],uuu,gg(i));
            gmax(i)=fabs(pfmc1[i]->gmax);
            if (!initrun_flag)
            {
              if (gmax(i)<1.e-4  || pfmc1[i]->ireturn<=0)
              {
                icon(i)=1;
              }
              else
              {
                converged=0;
              }
            }
          }
        }
      }
      initrun_flag=0;
      for (int i2=1;i2<=num_separable_calls;i2++)
      {
        int m=(*derindex)(i2).indexmax();
        for (int j=1;j<=m;j++)
        {
          u((*derindex)(i2)(j))=uu(i2,j);
        }
      }
      // put the
      //if (fmc1.ireturn>0)
      {
        dvariable vf=0.0;
        pen=initial_params::reset(dvar_vector(u));
        *objective_function_value::pobjfun=0.0;

        //num_separable_calls=0;

        pmin->inner_opt_flag=1;
        pfmin->AD_uf_inner();
        pmin->inner_opt_flag=0;

        if (saddlepointflag)
        {
          *objective_function_value::pobjfun*=-1.0;
        }
        if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
        {
          quadratic_prior::get_M_calculations();
        }
        vf+=*objective_function_value::pobjfun;

        objective_function_value::fun_without_pen=value(vf);
        vf+=pen;

        gradcalc(usize,g);
        dvector* pggi = &gg(1);
        for (int i=1;i<=num_separable_calls;i++)
        {
          int m=(*derindex)(i).indexmax();
          double* pggij = pggi->get_v() + 1;
          for (int j=1;j<=m;j++)
          {
            *pggij = g((*derindex)(i)(j));

            ++pggij;
          }

          ++pggi;
        }
  /*
        {
          ofstream ofs("l:/temp1.dat");
          ofs << g.indexmax() << " " << setprecision(15) << g << endl;
        }
  */
        if (saddlepointflag==2)
        {
          ff[1]=-(*separable_function_difference)(1);
          for (int i=2;i<=num_separable_calls;i++)
          {
            ff[i]=-(*separable_function_difference)(i);
            //ff[i]=-(*separable_function_difference)(i)
             // +(*separable_function_difference)(i-1);

            if (ff[i] < ffb[i])
            {
              ffb[i]=ff[i];
              uub[i]=uu[i];
              ggb[i]=gg[i];
            }
          }
        }
        else
        {
          ff[1]=(*separable_function_difference)(1);
          for (int i=2;i<=num_separable_calls;i++)
          {
            ff[i]=(*separable_function_difference)(i);
            //ff[i]=(*separable_function_difference)(i)
             // -(*separable_function_difference)(i-1);

            if (ff[i] < ffb[i])
            {
              ffb[i]=ff[i];
              uub[i]=uu[i];
              ggb[i]=gg[i];
            }
          }
        }
        f=0.0;
        for (int i2=1;i2<=num_separable_calls;i2++)
        {
          f+=ff[i2];
        }
        if (f<fb)
        {
          fb=f;
          ub=u;
        }
      }
      u=ub;
    }
    double tmax=max(gmax);
    output_stream <<  " inner maxg = " << std::scientific << setprecision(10) <<  tmax;
    if (tmax< 1.e-4) break;
  }
  fmc1.ireturn=0;
  fmc1.fbest=fb;
  gradient_structure::set_NO_DERIVATIVES();
  //num_separable_calls=0;
  pmin->inner_opt_flag=1;
  pfmin->AD_uf_inner();
  pmin->inner_opt_flag=0;
  if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
  {
    quadratic_prior::get_M_calculations();
  }
  gradient_structure::set_YES_DERIVATIVES();
  for (int i=1;i<=num_separable_calls;i++)
  {
    if (pfmc1[i])
    {
      delete pfmc1[i];
    }
  }
  pfmc1++;
  delete [] pfmc1;
  pfmc1 = 0;
  return u;
}