ADMB Documentation  11.1.2532
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_ma6.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fvar_ma6.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00013 void dfempirical_covarv_partial(void);
00014 
00019 dvar_matrix empirical_covariance(const dvar_matrix& _v1,
00020   const imatrix& _missflags)
00021  {
00022    dvar_matrix& v1 = (dvar_matrix&) (_v1);
00023    imatrix& missflags=(imatrix&) (_missflags);
00024    int mmin=v1(v1.indexmin()).indexmin();
00025    int mmax=v1(v1.indexmin()).indexmax();
00026    dvar_matrix tmp(mmin,mmax,mmin,mmax);
00027    int rmin=v1.indexmin();
00028    int rmax=v1.indexmax();
00029    int nobs=rmax-rmin+1;
00030 
00031    tmp.initialize();
00032    for (int ii=rmin; ii<=rmax; ii++)
00033    {
00034      for (int i=mmin; i<=mmax; i++)
00035      {
00036        if (!missflags(ii,i))
00037        {
00038          for (int j=mmin; j<=mmax; j++)
00039          {
00040            if (!missflags(ii,j))
00041            {
00042              tmp.elem_value(i,j)+=v1.elem_value(ii,i)*v1.elem_value(ii,j);
00043            }
00044          }
00045        }
00046      }
00047    }
00048    for (int i=mmin; i<=mmax; i++)
00049    {
00050      for (int j=mmin; j<=mmax; j++)
00051      {
00052        tmp.elem_value(i,j)/=nobs;
00053      }
00054    }
00055   save_identifier_string("ru");
00056   missflags.save_imatrix_value();
00057   missflags.save_imatrix_position();
00058   save_int_value(nobs);
00059   tmp.save_dvar_matrix_position();
00060   v1.save_dvar_matrix_value();
00061   v1.save_dvar_matrix_position();
00062   save_identifier_string("rv");
00063   gradient_structure::GRAD_STACK1->
00064       set_gradient_stack(dfempirical_covarv_partial);
00065    return(tmp);
00066  }
00067 
00072 void dfempirical_covarv_partial(void)
00073 {
00074   verify_identifier_string("rv");
00075   dvar_matrix_position v1pos=restore_dvar_matrix_position();
00076   dmatrix v1=restore_dvar_matrix_value(v1pos);
00077   dvar_matrix_position tmppos=restore_dvar_matrix_position();
00078   dmatrix dftmp=restore_dvar_matrix_derivatives(tmppos);
00079   int nobs=restore_int_value();
00080   imatrix_position mfpos=restore_imatrix_position();
00081   imatrix missflags=restore_imatrix_value(mfpos);
00082   verify_identifier_string("ru");
00083   int mmin=v1(v1.indexmin()).indexmin();
00084   int mmax=v1(v1.indexmin()).indexmax();
00085   int rmin=v1.indexmin();
00086   int rmax=v1.indexmax();
00087 
00088   dmatrix dfv1(rmin,rmax,mmin,mmax);
00089   dfv1.initialize();
00090   for (int i=mmin; i<=mmax; i++)
00091   {
00092     for (int j=mmin; j<=mmax; j++)
00093     {
00094       //tmp.elem_value(i,j)/=nobs;
00095       dftmp(i,j)/=nobs;
00096     }
00097   }
00098   for (int ii=rmin; ii<=rmax; ii++)
00099   {
00100     for (int i=mmin; i<=mmax; i++)
00101     {
00102       if (!missflags(ii,i))
00103       {
00104         for (int j=mmin; j<=mmax; j++)
00105         {
00106           if (!missflags(ii,j))
00107           {
00108             //tmp.elem_value(i,j)+=v1.elem_value(ii,i)*v1.elem_value(ii,j);
00109             dfv1(ii,i)+=dftmp(i,j)*v1(ii,j);
00110             dfv1(ii,j)+=dftmp(i,j)*v1(ii,i);
00111           }
00112         }
00113       }
00114     }
00115   }
00116   dfv1.save_dmatrix_derivatives(v1pos);
00117 }
00118 
00119 void dfempirical_covarv(void);
00120 
00125 dvar_matrix empirical_covariance(const dvar_matrix& v1)
00126  {
00127    int mmin=v1(v1.indexmin()).indexmin();
00128    int mmax=v1(v1.indexmin()).indexmax();
00129    dvar_matrix tmp(mmin,mmax,mmin,mmax);
00130    int rmin=v1.indexmin();
00131    int rmax=v1.indexmax();
00132 
00133    tmp.initialize();
00134    int nobs=rmax-rmin+1;
00135    for (int ii=rmin; ii<=rmax; ii++)
00136    {
00137      for (int i=mmin; i<=mmax; i++)
00138      {
00139        for (int j=mmin; j<=mmax; j++)
00140        {
00141          tmp.elem_value(i,j)+=v1.elem_value(ii,i)*v1.elem_value(ii,j);
00142        }
00143      }
00144    }
00145    for (int i=mmin; i<=mmax; i++)
00146    {
00147      for (int j=mmin; j<=mmax; j++)
00148      {
00149        tmp.elem_value(i,j)/=double(nobs);
00150      }
00151    }
00152   save_identifier_string("ru");
00153   tmp.save_dvar_matrix_position();
00154   v1.save_dvar_matrix_value();
00155   v1.save_dvar_matrix_position();
00156   save_identifier_string("rv");
00157   gradient_structure::GRAD_STACK1->
00158       set_gradient_stack(dfempirical_covarv);
00159    return(tmp);
00160  }
00161 
00166 void dfempirical_covarv(void)
00167 {
00168   verify_identifier_string("rv");
00169   dvar_matrix_position v1pos=restore_dvar_matrix_position();
00170   dmatrix v1=restore_dvar_matrix_value(v1pos);
00171   dvar_matrix_position tmppos=restore_dvar_matrix_position();
00172   dmatrix dftmp=restore_dvar_matrix_derivatives(tmppos);
00173   verify_identifier_string("ru");
00174   int mmin=v1(v1.indexmin()).indexmin();
00175   int mmax=v1(v1.indexmin()).indexmax();
00176   int rmin=v1.indexmin();
00177   int rmax=v1.indexmax();
00178 
00179   dmatrix dfv1(rmin,rmax,mmin,mmax);
00180   dfv1.initialize();
00181   int nobs=rmax-rmin+1;
00182   for (int i=mmin; i<=mmax; i++)
00183   {
00184     for (int j=mmin; j<=mmax; j++)
00185     {
00186       //tmp.elem_value(i,j)/=double(nobs);
00187       dftmp(i,j)/=double(nobs);
00188     }
00189   }
00190   for (int ii=rmin; ii<=rmax; ii++)
00191   {
00192     for (int i=mmin; i<=mmax; i++)
00193     {
00194       for (int j=mmin; j<=mmax; j++)
00195       {
00196         //tmp.elem_value(i,j)+=v1.elem_value(ii,i)*v1.elem_value(ii,j);
00197         dfv1(ii,i)+=dftmp(i,j)*v1(ii,j);
00198         dfv1(ii,j)+=dftmp(i,j)*v1(ii,i);
00199       }
00200     }
00201   }
00202   dfv1.save_dmatrix_derivatives(v1pos);
00203 }
00204 
00205 void dfouter_prodvv(void);
00206 
00211 dvar_matrix outer_prod(const dvar_vector& v1, const dvar_vector& v2)
00212  {
00213    dvar_matrix tmp(v1.indexmin(),v1.indexmax(), v2.indexmin(), v2.indexmax() );
00214 
00215    for (int i=v1.indexmin(); i<=v1.indexmax(); i++)
00216    {
00217      for (int j=v2.indexmin(); j<=v2.indexmax(); j++)
00218      {
00219        tmp.elem_value(i,j)=v1.elem_value(i)*v2.elem_value(j);
00220      }
00221    }
00222   save_identifier_string("tu");
00223   tmp.save_dvar_matrix_position();
00224   v1.save_dvar_vector_value();
00225   v1.save_dvar_vector_position();
00226   v2.save_dvar_vector_value();
00227   v2.save_dvar_vector_position();
00228   save_identifier_string("tv");
00229   gradient_structure::GRAD_STACK1->
00230       set_gradient_stack(dfouter_prodvv);
00231    return(tmp);
00232  }
00233 
00238 void dfouter_prodvv(void)
00239 {
00240   verify_identifier_string("tv");
00241   dvar_vector_position v2pos=restore_dvar_vector_position();
00242   dvector v2=restore_dvar_vector_value(v2pos);
00243   dvar_vector_position v1pos=restore_dvar_vector_position();
00244   dvector v1=restore_dvar_vector_value(v1pos);
00245   dvar_matrix_position tmppos=restore_dvar_matrix_position();
00246   dmatrix dftmp=restore_dvar_matrix_derivatives(tmppos);
00247   verify_identifier_string("tu");
00248   dvector dfv1(v1pos.indexmin(),v1pos.indexmax());
00249   dvector dfv2(v2pos.indexmin(),v2pos.indexmax());
00250   dfv1.initialize();
00251   dfv2.initialize();
00252   for (int i=v1.indexmin(); i<=v1.indexmax(); i++)
00253   {
00254     for (int j=v2.indexmin(); j<=v2.indexmax(); j++)
00255     {
00256       //tmp.elem_value(i,j)=v1.elem_value(i)*v2.elem_value(j);
00257       dfv1(i)+=dftmp(i,j)*v2(j);
00258       dfv2(j)+=dftmp(i,j)*v1(i);
00259     }
00260   }
00261   dfv1.save_dvector_derivatives(v1pos);
00262   dfv2.save_dvector_derivatives(v2pos);
00263 }
00264 
00269 dvar_matrix outer_prod(const dvector& v1, const dvar_vector& v2)
00270  {
00271    RETURN_ARRAYS_INCREMENT();
00272 
00273    dvar_matrix tmp(v1.indexmin(),v1.indexmax(), v2.indexmin(), v2.indexmax() );
00274 
00275    for (int i=v1.indexmin(); i<=v1.indexmax(); i++)
00276    {
00277      for (int j=v2.indexmin(); j<=v2.indexmax(); j++)
00278      {
00279        tmp.elem(i,j)=v1.elem(i)*v2.elem(j);
00280      }
00281    }
00282    RETURN_ARRAYS_DECREMENT();
00283    return(tmp);
00284  }
00285 
00290 dvar_matrix outer_prod(const dvar_vector& v1, const dvector& v2)
00291  {
00292    RETURN_ARRAYS_INCREMENT();
00293 
00294    dvar_matrix tmp(v1.indexmin(),v1.indexmax(), v2.indexmin(), v2.indexmax() );
00295 
00296    for (int i=v1.indexmin(); i<=v1.indexmax(); i++)
00297    {
00298      for (int j=v2.indexmin(); j<=v2.indexmax(); j++)
00299      {
00300        tmp.elem(i,j)=v1.elem(i)*v2.elem(j);
00301      }
00302    }
00303    RETURN_ARRAYS_DECREMENT();
00304    return(tmp);
00305  }