From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 24278 invoked by alias); 18 Mar 2004 09:15:15 -0000 Mailing-List: contact gcc-help-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Archive: List-Post: List-Help: Sender: gcc-help-owner@gcc.gnu.org Received: (qmail 24178 invoked from network); 18 Mar 2004 09:15:09 -0000 Received: from unknown (HELO filabres.ual.es) (150.214.156.2) by sources.redhat.com with SMTP; 18 Mar 2004 09:15:09 -0000 Received: from filabres.ual.es (hall.ual.es [150.214.156.3]) by filabres.ual.es (8.12.9/8.9.3) with ESMTP id i2I9F5hS003287 for ; Thu, 18 Mar 2004 10:15:06 +0100 (MET) Received: from dali.ace.ual.es (dali.ace.ual.es [193.147.118.56]) by filabres.ual.es (8.12.9/8.9.3) with ESMTP id i2I9F4DD003277; Thu, 18 Mar 2004 10:15:04 +0100 (MET) Received: from [193.147.118.93] ( [193.147.118.93]) by dali.ace.ual.es (8.11.6/8.11.6) with ESMTP id i2IAtlL09740; Thu, 18 Mar 2004 11:55:47 +0100 Subject: Re: Please, help with a weir error on templates. From: Jose Roman Bilbao To: "Alex J. Dam" Cc: gcc-help@gcc.gnu.org In-Reply-To: <200403180007.51842.alexjdam@yahoo.com.br> References: <1079534582.10524.9.camel@gris.ace.ual.es> <200403180007.51842.alexjdam@yahoo.com.br> Content-Type: multipart/mixed; boundary="=-vj1wM3x7pj21v3gE0Vq3" Message-Id: <1079601306.13642.6.camel@gris.ace.ual.es> Mime-Version: 1.0 Date: Thu, 18 Mar 2004 15:19:00 -0000 X-SW-Source: 2004-03/txt/msg00156.txt.bz2 --=-vj1wM3x7pj21v3gE0Vq3 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Content-length: 1166 Hi, As you asked I attatch here the two files which generate the problem. One is the .inc file where code is implemented and .cc and .hh files which make use of the .inc file. I am using gcc coming with Mandrake 10 3.3.2 and this is the g++ call: g++ -DHAVE_CONFIG_H -I. -I. -I../.. -I../../Lib -I../../Lib/XmippData/Bilib/headers -I../../Lib/XmippData/Bilib/types -I../../Lib/XmippData/Bilib -g -O2 -MT xmippMatrices1D.lo -MD -MP -MF .deps/xmippMatrices1D.Tpo -c Src/xmippMatrices1D.cc -fPIC -DPIC -o .libs/xmippMatrices1D.lo Thanks El jue, 18-03-2004 a las 04:07, Alex J. Dam escribió: > Can you post a whole C++ file which generates that error? > > Also, which GCC version are you using? Which options are you passing to g++? > > > Jose Roman Bilbao scripsit: > > When compiler reaches the next code: > > > > template <> > > void core_array_by_scalar< complex >(const maTC &op1, > > const complex &op2, maTC &result, char operation) _THROW { > > > > it complains with: > > > > Src/MultidimBasic.inc:332: error: `double_complex' was not declared in > > this scope > > Src/MultidimBasic.inc:332: error: template argument 1 is invalid --=-vj1wM3x7pj21v3gE0Vq3 Content-Disposition: attachment; filename=MultidimBasic.inc Content-Type: text/plain; name=MultidimBasic.inc; charset=UTF-8 Content-Transfer-Encoding: 7bit Content-length: 16715 /*************************************************************************** * * Authors: Carlos Oscar S. Sorzano (coss@cnb.uam.es) * * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC * * 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 2 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 Aouble_ 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, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA * 02111-1307 USA * * All comments concerning this program package may be sent to the * e-mail address 'xmipp@cnb.uam.es' ***************************************************************************/ /* ------------------------------------------------------------------------- */ /* MULTIDIM BASIC */ /* ------------------------------------------------------------------------- */ #include #define mi MULTIDIM_ELEM(*this,i) #define msize MULTIDIM_SIZE(*this) /* Print stats ------------------------------------------------------------- */ template void maT::print_stats(ostream &out) const { T min_val, max_val; double avg_val, dev_val; compute_stats(avg_val, dev_val, min_val, max_val); out.setf(ios::showpoint); int old_prec=out.precision(7); out << " min= "; out.width(9); out << min_val; out << " max= "; out.width(9); out << max_val; out << " avg= "; out.width(9); out << avg_val; out << " dev= "; out.width(9); out << dev_val; out.precision(old_prec); } /* Compute max ------------------------------------------------------------- */ template T maT::compute_max() const { if (__dim<=0) return (T)0; T max_val=MULTIDIM_ELEM(*this,0); FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) if (mi>max_val) max_val=mi; return max_val; } /* Compute min ------------------------------------------------------------- */ template T maT::compute_min() const { if (__dim<=0) return (T)0; T min_val=MULTIDIM_ELEM(*this,0); FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) if (mi void maT::compute_double_minmax(double &min, double &max) const { min=max=0.0f; if (__dim<=0) return; min=max=(double) MULTIDIM_ELEM(*this,0); FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) { if (mimax) max=(double)mi; } } template void maT::compute_double_minmax(double &min, double &max, const matrix1D &corner1, const matrix1D &corner2) const { matrix1D dcorner1, dcorner2; type_cast(corner1,dcorner1); type_cast(corner2,dcorner2); compute_double_minmax(min,max,dcorner1,dcorner2); } template void maT::compute_stats(double &avg, double &stddev, T &min_val, T &max_val, const matrix1D &corner1, const matrix1D &corner2) const { matrix1D dcorner1, dcorner2; type_cast(corner1,dcorner1); type_cast(corner2,dcorner2); compute_stats(avg,stddev,min_val,max_val,dcorner1,dcorner2); } /* Compute avg ------------------------------------------------------------- */ template double maT::compute_avg() const { if (__dim<=0) return 0; double sum=0; FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) sum += (double) mi; return sum/msize; } /* Compute stddev ---------------------------------------------------------- */ template double maT::compute_stddev() const { if (msize<=0) return 0; double avg=0, stddev=0; FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) { avg += (double) mi; stddev += (double) mi * (double) mi; } if (msize>1) { avg/=msize; stddev = stddev/msize - avg*avg; stddev*= msize/(msize-1); stddev = sqrt((double)(ABS(stddev))); // Foreseeing numerical // instabilities } else stddev=0; return stddev; } /* Compute stats ---------------------------------------------------------- */ template void maT::compute_stats(double &avg, double &stddev, T &min, T &max) const { avg=0; stddev=0; min=(T)0; max=(T)0; if (msize<=0) return; min=max=MULTIDIM_ELEM(*this,0); FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) { avg += (double) mi; stddev += (double) mi * (double) mi ; if (mi>max) max=mi; if (mi1) { stddev = stddev/msize - avg*avg; stddev*= msize/(msize-1); stddev = sqrt((double)(ABS(stddev))); // Foreseeing numerical // instabilities } else stddev=0; } /* Range adjust ------------------------------------------------------------ */ // This function takes a vector with range between min0 ... max0 and // linearly transforms it to minF ... maxF. // If the input vector is a constant then it is adjusted to minF template void maT::range_adjust(T minF, T maxF) { if (msize==0) return; T min0=compute_min(); T max0=compute_max(); // If max0==min0, it means that the vector is a constant one, so the // only possible transformation is to a fixed minF double slope; if (max0!=min0) slope=(double)(maxF-minF)/(double)(max0-min0); else slope=0; FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) mi=minF+(T) (slope * (double) (mi-min0)); } /* Statistics Adjust ------------------------------------------------------- */ // This function takes a vector with an average avg0 and standard deviation // dev0, and transforms it linearly into a vector with avgF and devF // If input vector has got dev0==0, ie, is a constant vector then only the // average is adjusted, and the deviation remains being 0. template void maT::statistics_adjust(double avgF, double stddevF) { double avg0,stddev0; double a,b; if (msize==0) return; T min, max; compute_stats(avg0, stddev0, min, max); if (stddev0!=0) a=(double)stddevF/(double)stddev0; else a=0; b=(double)avgF - a*(double)avg0; FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) mi = (T) (a*(double)mi+b); } /* Effective range --------------------------------------------------------- */ template double maT::effective_range(double percentil_out) { // histogram1D hist; // compute_hist(*this,hist,200); // double min_val = hist.percentil(percentil_out/2); // double max_val = hist.percentil(100-percentil_out/2); // return max_val-min_val; return 0; } /* Init random ------------------------------------------------------------- */ template void maT::init_random(double op1, double op2, const string &mode) _THROW { if (mode=="uniform") FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) mi=(T) rnd_unif(op1,op2); else if (mode=="gaussian") FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) mi=(T) rnd_gaus(op1,op2); else REPORT_ERROR(1005,(string)"Init_random: Mode not supported ("+mode+")"); } /* Add Noise --------------------------------------------------------------- */ // Supported random distributions: // "uniform", between op1 and op2 // "gaussian", with avg=op1 and var=op2 // It is not an error that the vector is empty template void maT::add_noise(double op1, double op2, const string &mode) const _THROW { if (mode=="uniform") FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) mi +=(T) rnd_unif(op1,op2); else if (mode=="gaussian") FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) mi +=(T) rnd_gaus(op1,op2); else REPORT_ERROR(1005,(string)"Add_noise: Mode not supported ("+mode+")"); } /* Threshold --------------------------------------------------------------- */ // Substitute component values by other according to the type of threshold // to apply // abs_above, abs_below, above, below, range template void maT::threshold(const string &type, T a, T b) { int mode; if (type == "abs_above") mode=1; else if (type == "abs_below") mode=2; else if (type == "above") mode=3; else if (type == "below") mode=4; else if (type == "range") mode=5; else REPORT_ERROR(1005, (string)"Threshold: mode not supported ("+type+")"); FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) switch (mode) { case 1: if (ABS(mi)>a) mi=SGN(mi)*b; break; case 2: if (ABS(mi)a) mi=b; break; case 4: if (mib) mi=b; break; } } /* Count with threshold ---------------------------------------------------- */ template long maT::count_threshold(const string &type, T a, T b) { int mode; if (type == "abs_above") mode=1; else if (type == "abs_below") mode=2; else if (type == "above") mode=3; else if (type == "below") mode=4; else if (type == "range") mode=5; else REPORT_ERROR(1005, (string)"Threshold: mode not supported ("+type+")"); long retval=0; FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(*this) switch (mode) { case 1: if (ABS(mi)>a) retval++; break; case 2: if (ABS(mi)a) retval++; break; case 4: if (mi=a && mi<=b) retval++; break; } return retval; } /* Equality for normal data types ------------------------------------------ */ template bool operator == (const maT &op1, const maT &op2) {return op1.equal(op2);} /* Equality for complex numbers -------------------------------------------- */ #ifdef ENABLE_COMPLEX_MATRICES bool operator == (const ma< complex > &op1, const ma< complex > &op2) { if (!op1.same_shape(op2)) return FALSE; FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(op1) if (ABS(MULTIDIM_ELEM(op1,i).real()-MULTIDIM_ELEM(op2,i).real()) >XMIPP_EQUAL_ACCURACY ||ABS(MULTIDIM_ELEM(op1,i).imag()-MULTIDIM_ELEM(op2,i).imag()) >XMIPP_EQUAL_ACCURACY) return FALSE; return TRUE; } #endif /* Core array by scalar ---------------------------------------------------- */ template void core_array_by_scalar(const maT &op1, const T &op2, maT &result, char operation) _THROW { FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(result) switch (operation) { case '+': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) + op2; break; case '-': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) - op2; break; case '*': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) * op2; break; case '/': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) / op2; break; case '^': MULTIDIM_ELEM(result,i)= (T) pow((double)MULTIDIM_ELEM(op1,i),(double)op2); break; } } template <> void core_array_by_scalar< complex >(const maTC &op1, const complex &op2, maTC &result, char operation) _THROW { FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(result) switch (operation) { case '+': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) + op2; break; case '-': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) - op2; break; case '*': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) * op2; break; case '/': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) / op2; break; case '^': MULTIDIM_ELEM(result,i)= pow(MULTIDIM_ELEM(op1,i),op2); break; } } /* Scalar by array --------------------------------------------------------- */ template void core_scalar_by_array(const T &op1, const maT &op2, maT &result, char operation) _THROW { FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(result) switch (operation) { case '+': MULTIDIM_ELEM(result,i)=op1 + MULTIDIM_ELEM(op2,i); break; case '-': MULTIDIM_ELEM(result,i)=op1 - MULTIDIM_ELEM(op2,i); break; case '*': MULTIDIM_ELEM(result,i)=op1 * MULTIDIM_ELEM(op2,i); break; case '/': MULTIDIM_ELEM(result,i)=op1 / MULTIDIM_ELEM(op2,i); break; case '^': MULTIDIM_ELEM(result,i)=(T) pow((double)op1,(double)MULTIDIM_ELEM(op2,i)); break; } } template <> void core_scalar_by_array< complex >(const complex &op1, const maTC &op2, maTC &result, char operation) _THROW { FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(result) switch (operation) { case '+': MULTIDIM_ELEM(result,i)=op1 + MULTIDIM_ELEM(op2,i); break; case '-': MULTIDIM_ELEM(result,i)=op1 - MULTIDIM_ELEM(op2,i); break; case '*': MULTIDIM_ELEM(result,i)=op1 * MULTIDIM_ELEM(op2,i); break; case '/': MULTIDIM_ELEM(result,i)=op1 / MULTIDIM_ELEM(op2,i); break; case '^': MULTIDIM_ELEM(result,i)=pow(op1,MULTIDIM_ELEM(op2,i)); break; } } /* Array by array ---------------------------------------------------------- */ template void core_array_by_array(const maT &op1, const maT &op2, maT &result, char operation) _THROW { FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(result) switch (operation) { case '+': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) + MULTIDIM_ELEM(op2,i); break; case '-': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) - MULTIDIM_ELEM(op2,i); break; case '*': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) * MULTIDIM_ELEM(op2,i); break; case '/': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) / MULTIDIM_ELEM(op2,i); break; case '^': MULTIDIM_ELEM(result,i)= (T) pow((double)MULTIDIM_ELEM(op1,i), (double)MULTIDIM_ELEM(op2,i)); break; } } template <> void core_array_by_array< complex >(const maTC &op1, const maTC &op2, maTC &result, char operation) _THROW { FOR_ALL_ELEMENTS_IN_MULTIDIM_ARRAY(result) switch (operation) { case '+': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) + MULTIDIM_ELEM(op2,i); break; case '-': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) - MULTIDIM_ELEM(op2,i); break; case '*': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) * MULTIDIM_ELEM(op2,i); break; case '/': MULTIDIM_ELEM(result,i)=MULTIDIM_ELEM(op1,i) / MULTIDIM_ELEM(op2,i); break; case '^': MULTIDIM_ELEM(result,i)=pow(MULTIDIM_ELEM(op1,i), MULTIDIM_ELEM(op2,i)); break; } } --=-vj1wM3x7pj21v3gE0Vq3 Content-Disposition: attachment; filename=xmippMatrices1D.cc Content-Type: text/x-c++; name=xmippMatrices1D.cc; charset=UTF-8 Content-Transfer-Encoding: 7bit Content-length: 17192 /*************************************************************************** * * Authors: Carlos Oscar S. Sorzano (coss@cnb.uam.es) * * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC * * 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 2 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, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA * 02111-1307 USA * * All comments concerning this program package may be sent to the * e-mail address 'xmipp@cnb.uam.es' ***************************************************************************/ /* ------------------------------------------------------------------------- */ /* VECTORS */ /* ------------------------------------------------------------------------- */ #include "../xmippMatrices1D.hh" #include "../xmippArgs.hh" #include /* ************************************************************************* */ /* IMPLEMENTATIONS */ /* ************************************************************************* */ #define maT matrix1D #define ma matrix1D #include "MultidimBasic.inc" #undef ma #undef maT // This file contain several general functions for arithmetic operations // They are defined for a general multidimensional_array, but // we have "redirected" maT to matrix1D, this way they will do OK in // this vector library /* Print shape ------------------------------------------------------------- */ template void vT::print_shape(ostream &out) const { out << "Size: " << XSIZE(*this) << "i=[" << STARTINGX(*this) << ".." << FINISHINGX(*this) << "]"; } /* Get size--- ------------------------------------------------------------- */ template void vT::get_size(int *size) const {size[0]=xdim; size[1]=1; size[2]=1;} /* Outside ----------------------------------------------------------------- */ template bool vT::outside(const matrix1D &v) const _THROW { if (XSIZE(v)<1) REPORT_ERROR(1,"Outside: index vector has got not enough components"); return (XX(v)FINISHINGX(*this)); } template bool vT::outside(int i) const { return (iFINISHINGX(*this)); } /* Intersects -------------------------------------------------------------- */ template bool vT::intersects(const vT &m) const {return intersects(STARTINGX(m), XSIZE(m)-1);} template bool vT::intersects(const matrix1D &corner1, const matrix1D &corner2) const _THROW { if (XSIZE(corner1)!=1 || XSIZE(corner2)!=1) REPORT_ERROR(1002,"intersects 1D: corner sizes are not 1"); return intersects(XX(corner1),XX(corner2)-XX(corner1)); } template bool vT::intersects(double x0, double xdim) const { SPEED_UP_temps; spduptmp0=MAX(STARTINGX(*this), x0); spduptmp1=MIN(FINISHINGX(*this),x0+xdim); if (spduptmp0>spduptmp1) return FALSE; return TRUE; } /* Corner ------------------------------------------------------------------ */ template bool vT::isCorner(const matrix1D &v) _THROW { if (XSIZE(v)<1) REPORT_ERROR(1,"isCorner: index vector has got not enough components"); return (XX(v)==STARTINGX(*this) || XX(v)==FINISHINGX(*this)); } /* Border ------------------------------------------------------------------ */ template bool vT::isBorder(const matrix1D &v) _THROW { if (XSIZE(v)<1) REPORT_ERROR(1,"isBorder: index vector has got not enough components"); return isBorder(XX(v)); } template bool vT::isBorder(int i) { return (i==STARTINGX(*this) || i==FINISHINGX(*this)); } /* Patch ------------------------------------------------------------------- */ template void vT::patch(const vT &patch_array, char operation) { SPEED_UP_temps; FOR_ALL_ELEMENTS_IN_COMMON_IN_MATRIX1D(patch_array,*this) switch (operation) { case '=': VEC_ELEM(*this,i) =VEC_ELEM(patch_array,i); break; case '+': VEC_ELEM(*this,i)+=VEC_ELEM(patch_array,i); break; case '-': VEC_ELEM(*this,i)-=VEC_ELEM(patch_array,i); break; case '*': VEC_ELEM(*this,i)*=VEC_ELEM(patch_array,i); break; case '/': VEC_ELEM(*this,i)/=VEC_ELEM(patch_array,i); break; } } /* Output stream ----------------------------------------------------------- */ template ostream& operator << (ostream& out, const vT& v) { if (MULTIDIM_SIZE(v)==0) out << "NULL vector\n"; else { // Look for the exponent vT aux(v); aux.ABSnD(); int prec=best_prec(aux.compute_max(),10); FOR_ALL_ELEMENTS_IN_MATRIX1D(v) { if (v.row) out << FtoA((double)VEC_ELEM(v,i),10,prec) << ' '; else out << FtoA((double)VEC_ELEM(v,i),10,prec) << '\n'; } } return out; } // Special case for complex numbers template <> ostream& operator << (ostream& out, const matrix1D< complex >& v) { if (MULTIDIM_SIZE(v)==0) out << "NULL vector\n"; else { FOR_ALL_ELEMENTS_IN_MATRIX1D(v) { if (v.row) out << VEC_ELEM(v,i) << ' '; else out << VEC_ELEM(v,i) << '\n'; } } return out; } /* Linear initialisation --------------------------------------------------- */ // It is not an error that the vector is empty template void vT::init_linear(T minF, T maxF, int n, const string &mode) _THROW { double slope; int steps; if (mode=="incr") { steps=1+(int) FLOOR((double)ABS((maxF-minF))/((double) n)); slope=n*SGN(maxF-minF); } else if (mode=="steps") { steps=n; slope=(maxF-minF)/(steps-1); } else REPORT_ERROR(1005,"Init_linear: Mode not supported ("+mode+")"); if (steps==0) clear(); else { resize(steps); for (int i=0; i void vT::self_reverse() { int imax=(int)(XSIZE(*this)-1)/2; for (int i=0; i<=imax; i++) { T aux; SWAP(MULTIDIM_ELEM(*this,i),MULTIDIM_ELEM(*this,XSIZE(*this)-1-i),aux); } } /* Sort vector ------------------------------------------------------------- */ template vT vT::sort() const { vT temp; matrix1D aux; if (XSIZE(*this)==0) return temp; // Initialise data type_cast(*this, aux); // Sort double * aux_array=aux.adapt_for_numerical_recipes(); qcksrt(xdim,aux_array); type_cast(aux,temp); return temp; } /* Index in order ---------------------------------------------------------- */ template matrix1D vT::index_sort() const { matrix1D indx; matrix1D temp; if (XSIZE(*this)==0) return indx; if (XSIZE(*this)==1) { indx.resize(1); indx(0)=1; return indx; } // Initialise data indx.resize(xdim); type_cast(*this,temp); // Sort indexes double * temp_array=temp.adapt_for_numerical_recipes(); int * indx_array=indx.adapt_for_numerical_recipes(); indexx(xdim,temp_array,indx_array); return indx; } /* Window ------------------------------------------------------------------ */ template void vT::window(int x0, int xF, T init_value) { vT result(xF-x0+1); STARTINGX(result)=x0; for (int j=x0; j<=xF; j++) if (j>=STARTINGX(*this) && j<=FINISHINGX(*this)) VEC_ELEM(result,j)=VEC_ELEM(*this,j); else VEC_ELEM(result,j)=init_value; *this=result; } /* Max index --------------------------------------------------------------- */ template void vT::max_index(int &imax) const { if (XSIZE(*this)==0) {imax=-1; return;} imax=0; T max=VEC_ELEM(*this,imax); FOR_ALL_ELEMENTS_IN_MATRIX1D(*this) if (VEC_ELEM(*this,i)>max) {max=VEC_ELEM(*this,i); imax=i;} } /* Min index --------------------------------------------------------------- */ template void vT::min_index(int &imin) const { if (XSIZE(*this)==0) {imin=-1; return;} imin=0; T min=VEC_ELEM(*this,imin); FOR_ALL_ELEMENTS_IN_MATRIX1D(*this) if (VEC_ELEM(*this,i) void vT::compute_stats(double &avg, double &stddev, T &min_val, T &max_val, const matrix1D &corner1, const matrix1D &corner2) const { min_val=max_val=(*this)(corner1); matrix1D r(3); double N=0, sum=0, sum2=0; FOR_ALL_ELEMENTS_IN_MATRIX1D_BETWEEN(corner1,corner2) { sum+=(*this)(r); sum2+=(*this)(r)*(*this)(r); N++; if ((*this)(r)max_val) max_val=(*this)(r); } if (N!=0) { avg=sum/N; stddev=sqrt(sum2/N-avg*avg); } else {avg=stddev=0;} } /* Minimum and maximum in region ------------------------------------------- */ template void vT::compute_double_minmax(double &min_val, double &max_val, const matrix1D &corner1, const matrix1D &corner2) const { min_val=max_val=(*this)(corner1); matrix1D r(1); FOR_ALL_ELEMENTS_IN_MATRIX1D_BETWEEN(corner1,corner2) { if ((*this)(r)max_val) max_val=(*this)(r); } } /* Vector R2 and R3 -------------------------------------------------------- */ matrix1D vector_R2(double x, double y) { matrix1D result(2); VEC_ELEM(result,0)=x; VEC_ELEM(result,1)=y; return result; } matrix1D vector_R3(double x, double y, double z) { matrix1D result(3); VEC_ELEM(result,0)=x; VEC_ELEM(result,1)=y; VEC_ELEM(result,2)=z; return result; } matrix1D vector_R3(int x, int y, int z) { matrix1D result(3); VEC_ELEM(result,0)=x; VEC_ELEM(result,1)=y; VEC_ELEM(result,2)=z; return result; } /* Are orthogonal ---------------------------------------------------------- */ int are_orthogonal (matrix1D &v1, matrix1D &v2, matrix1D &v3) _THROW { if (XSIZE(v1)!=3 || XSIZE(v2)!=3 || XSIZE(v3)!=3) REPORT_ERROR(1002,"Are orthogonal: Some vector do not belong to R3"); try{ if (dot_product(v1,v2)!=0) return 0; if (dot_product(v2,v3)!=0) return 0; if (dot_product(v1,v3)!=0) return 0; } catch (Xmipp_error) { REPORT_ERROR(1007,"Are orthogonal: Vectors are not all of the same shape"); } return 1; } /* Are system? ------------------------------------------------------------- */ int are_system (matrix1D &v1, matrix1D &v2, matrix1D &v3) _THROW { matrix1D aux(3); if (XSIZE(v1)!=3 || XSIZE(v2)!=3 || XSIZE(v3)!=3) REPORT_ERROR(1002,"Are orthogonal: Some vector do not belong to R3"); aux=vector_product(v1,v2); if (aux!=v3) return 0; aux=vector_product(v2,v3); if (aux!=v1) return 0; aux=vector_product(v3,v1); if (aux!=v2) return 0; return 1; } /* Cut to common size ------------------------------------------------------ */ template void cut_to_common_size(vT &V1, vT &V2) { int x0=MAX(STARTINGX(V1) ,STARTINGX(V2)); int xF=MIN(FINISHINGX(V1),FINISHINGX(V2)); V1.window(x0,xF); V2.window(x0,xF); } /* Sort two vectors -------------------------------------------------------- */ template void sort_two_vectors(vT &v1, vT &v2) _THROW { T temp; if (XSIZE(v1)!=XSIZE(v2) || STARTINGX(v1)!=STARTINGX(v2)) REPORT_ERROR(1007, "sort_two_vectors: vectors are not of the same shape"); FOR_ALL_ELEMENTS_IN_MATRIX1D(v1) { temp=MIN(VEC_ELEM(v1,i),VEC_ELEM(v2,i)); VEC_ELEM(v2,i)=MAX(VEC_ELEM(v1,i),VEC_ELEM(v2,i)); VEC_ELEM(v1,i)=temp; } } /* Powell's optimizer ------------------------------------------------------ */ void Powell_optimizer(matrix1D &p, int i0, int n, double (*f)(double *x), double ftol, double &fret, int &iter, const matrix1D &steps, bool show) { double *xi=NULL; // Adapt indexes of p double *pptr=p.adapt_for_numerical_recipes(); double *auxpptr=pptr+(i0-1); // Form direction matrix ask_Tvector(xi,1,n*n); int ptr; for (int i=1,ptr=1; i<=n; i++) for (int j=1; j<=n; j++,ptr++) xi[ptr]=(i==j)?steps(i-1):0; // Optimize xi-=n; // This is because NR works with matrices // starting at [1,1] powell(auxpptr,xi,n,ftol,iter,fret,f, show); xi+=n; // Exit free_Tvector(xi,1,n*n); } /* Center of mass ---------------------------------------------------------- */ template void vT::center_of_mass(matrix1D ¢er, void * mask) { center.init_zeros(1); double mass=0; matrix1D *imask=(matrix1D *) mask; FOR_ALL_ELEMENTS_IN_MATRIX1D(*this) { if (imask==NULL || VEC_ELEM(*imask,i)) { XX(center)+=i*VEC_ELEM(*this,i); mass+=VEC_ELEM(*this,i); } } if (mass!=0) center/=mass; } /* The multidimensional arrays must be instantiated if we want them to compile into a library. All methods of all possible classes must be called. The only basic types considered are: short int double double complex For complex arrays not all methods are valid, so they are not in the library. */ /* ------------------------------------------------------------------------- */ /* INSTANTIATE */ /* ------------------------------------------------------------------------- */ template void instantiate_vector (matrix1D v) { matrix1D a; matrix1D r; matrix1D ir; double d; T Taux; // General functions for multidimensional arrays a==a; a!=a; a=1-a; a=a-1; a=a*a; a.print_stats(); a.compute_double_minmax(d,d); a.compute_double_minmax(d,d,r,r); a.compute_double_minmax(d,d,ir,ir); a.compute_stats(d,d,Taux,Taux,ir,ir); a.compute_max(); a.compute_min(); a.compute_avg(); a.compute_stddev(); a.range_adjust(0,1); a.statistics_adjust(0,1); a.effective_range(99); a.init_random(0,1); a.add_noise(0,1); a.threshold("abs_above",1,0); a.count_threshold("abs_above",1,0); a.center_of_mass(r); a.print_shape(); int size[3]; a.get_size(size); a.outside(r); a.isBorder(0); matrix1D pixel(1); a.isBorder(pixel); a.outside(0); a.intersects(a); a.intersects(r,r); a.isCorner(r); a.patch(a); cout << a; a.window(0,1); int imax; a.max_index(imax); a.min_index(imax); cut_to_common_size(a,a); // Specific for vectors matrix1D indx; a.init_linear(0,5); a.reverse(); a=a.sort(); indx=a.index_sort(); sort_two_vectors(a,a); } void instantiate_complex_vector () { matrix1D< complex > a; matrix1D r; // General functions for multidimensional arrays a=(complex)1.0-a; a=a-(complex)1.0; a=a*a; a.print_shape(); int size[3]; a.get_size(size); a.outside(r); a.outside(0); a.intersects(a); a.intersects(r,r); a.isCorner(r); a.patch(a); cout << a; a.window(0,1); cut_to_common_size(a,a); // Specific for vectors a.reverse(); } void instantiate1D() { matrix1D V0; instantiate_vector(V0); matrix1D V1; instantiate_vector(V1); matrix1D V2; instantiate_vector(V2); matrix1D V3; instantiate_vector(V3); matrix1D V4; instantiate_vector(V4); matrix1D V5; instantiate_vector(V5); } --=-vj1wM3x7pj21v3gE0Vq3 Content-Disposition: attachment; filename=xmippMatrices1D.cc Content-Type: text/x-c++; name=xmippMatrices1D.cc; charset=UTF-8 Content-Transfer-Encoding: 7bit Content-length: 17192 /*************************************************************************** * * Authors: Carlos Oscar S. Sorzano (coss@cnb.uam.es) * * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC * * 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 2 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, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA * 02111-1307 USA * * All comments concerning this program package may be sent to the * e-mail address 'xmipp@cnb.uam.es' ***************************************************************************/ /* ------------------------------------------------------------------------- */ /* VECTORS */ /* ------------------------------------------------------------------------- */ #include "../xmippMatrices1D.hh" #include "../xmippArgs.hh" #include /* ************************************************************************* */ /* IMPLEMENTATIONS */ /* ************************************************************************* */ #define maT matrix1D #define ma matrix1D #include "MultidimBasic.inc" #undef ma #undef maT // This file contain several general functions for arithmetic operations // They are defined for a general multidimensional_array, but // we have "redirected" maT to matrix1D, this way they will do OK in // this vector library /* Print shape ------------------------------------------------------------- */ template void vT::print_shape(ostream &out) const { out << "Size: " << XSIZE(*this) << "i=[" << STARTINGX(*this) << ".." << FINISHINGX(*this) << "]"; } /* Get size--- ------------------------------------------------------------- */ template void vT::get_size(int *size) const {size[0]=xdim; size[1]=1; size[2]=1;} /* Outside ----------------------------------------------------------------- */ template bool vT::outside(const matrix1D &v) const _THROW { if (XSIZE(v)<1) REPORT_ERROR(1,"Outside: index vector has got not enough components"); return (XX(v)FINISHINGX(*this)); } template bool vT::outside(int i) const { return (iFINISHINGX(*this)); } /* Intersects -------------------------------------------------------------- */ template bool vT::intersects(const vT &m) const {return intersects(STARTINGX(m), XSIZE(m)-1);} template bool vT::intersects(const matrix1D &corner1, const matrix1D &corner2) const _THROW { if (XSIZE(corner1)!=1 || XSIZE(corner2)!=1) REPORT_ERROR(1002,"intersects 1D: corner sizes are not 1"); return intersects(XX(corner1),XX(corner2)-XX(corner1)); } template bool vT::intersects(double x0, double xdim) const { SPEED_UP_temps; spduptmp0=MAX(STARTINGX(*this), x0); spduptmp1=MIN(FINISHINGX(*this),x0+xdim); if (spduptmp0>spduptmp1) return FALSE; return TRUE; } /* Corner ------------------------------------------------------------------ */ template bool vT::isCorner(const matrix1D &v) _THROW { if (XSIZE(v)<1) REPORT_ERROR(1,"isCorner: index vector has got not enough components"); return (XX(v)==STARTINGX(*this) || XX(v)==FINISHINGX(*this)); } /* Border ------------------------------------------------------------------ */ template bool vT::isBorder(const matrix1D &v) _THROW { if (XSIZE(v)<1) REPORT_ERROR(1,"isBorder: index vector has got not enough components"); return isBorder(XX(v)); } template bool vT::isBorder(int i) { return (i==STARTINGX(*this) || i==FINISHINGX(*this)); } /* Patch ------------------------------------------------------------------- */ template void vT::patch(const vT &patch_array, char operation) { SPEED_UP_temps; FOR_ALL_ELEMENTS_IN_COMMON_IN_MATRIX1D(patch_array,*this) switch (operation) { case '=': VEC_ELEM(*this,i) =VEC_ELEM(patch_array,i); break; case '+': VEC_ELEM(*this,i)+=VEC_ELEM(patch_array,i); break; case '-': VEC_ELEM(*this,i)-=VEC_ELEM(patch_array,i); break; case '*': VEC_ELEM(*this,i)*=VEC_ELEM(patch_array,i); break; case '/': VEC_ELEM(*this,i)/=VEC_ELEM(patch_array,i); break; } } /* Output stream ----------------------------------------------------------- */ template ostream& operator << (ostream& out, const vT& v) { if (MULTIDIM_SIZE(v)==0) out << "NULL vector\n"; else { // Look for the exponent vT aux(v); aux.ABSnD(); int prec=best_prec(aux.compute_max(),10); FOR_ALL_ELEMENTS_IN_MATRIX1D(v) { if (v.row) out << FtoA((double)VEC_ELEM(v,i),10,prec) << ' '; else out << FtoA((double)VEC_ELEM(v,i),10,prec) << '\n'; } } return out; } // Special case for complex numbers template <> ostream& operator << (ostream& out, const matrix1D< complex >& v) { if (MULTIDIM_SIZE(v)==0) out << "NULL vector\n"; else { FOR_ALL_ELEMENTS_IN_MATRIX1D(v) { if (v.row) out << VEC_ELEM(v,i) << ' '; else out << VEC_ELEM(v,i) << '\n'; } } return out; } /* Linear initialisation --------------------------------------------------- */ // It is not an error that the vector is empty template void vT::init_linear(T minF, T maxF, int n, const string &mode) _THROW { double slope; int steps; if (mode=="incr") { steps=1+(int) FLOOR((double)ABS((maxF-minF))/((double) n)); slope=n*SGN(maxF-minF); } else if (mode=="steps") { steps=n; slope=(maxF-minF)/(steps-1); } else REPORT_ERROR(1005,"Init_linear: Mode not supported ("+mode+")"); if (steps==0) clear(); else { resize(steps); for (int i=0; i void vT::self_reverse() { int imax=(int)(XSIZE(*this)-1)/2; for (int i=0; i<=imax; i++) { T aux; SWAP(MULTIDIM_ELEM(*this,i),MULTIDIM_ELEM(*this,XSIZE(*this)-1-i),aux); } } /* Sort vector ------------------------------------------------------------- */ template vT vT::sort() const { vT temp; matrix1D aux; if (XSIZE(*this)==0) return temp; // Initialise data type_cast(*this, aux); // Sort double * aux_array=aux.adapt_for_numerical_recipes(); qcksrt(xdim,aux_array); type_cast(aux,temp); return temp; } /* Index in order ---------------------------------------------------------- */ template matrix1D vT::index_sort() const { matrix1D indx; matrix1D temp; if (XSIZE(*this)==0) return indx; if (XSIZE(*this)==1) { indx.resize(1); indx(0)=1; return indx; } // Initialise data indx.resize(xdim); type_cast(*this,temp); // Sort indexes double * temp_array=temp.adapt_for_numerical_recipes(); int * indx_array=indx.adapt_for_numerical_recipes(); indexx(xdim,temp_array,indx_array); return indx; } /* Window ------------------------------------------------------------------ */ template void vT::window(int x0, int xF, T init_value) { vT result(xF-x0+1); STARTINGX(result)=x0; for (int j=x0; j<=xF; j++) if (j>=STARTINGX(*this) && j<=FINISHINGX(*this)) VEC_ELEM(result,j)=VEC_ELEM(*this,j); else VEC_ELEM(result,j)=init_value; *this=result; } /* Max index --------------------------------------------------------------- */ template void vT::max_index(int &imax) const { if (XSIZE(*this)==0) {imax=-1; return;} imax=0; T max=VEC_ELEM(*this,imax); FOR_ALL_ELEMENTS_IN_MATRIX1D(*this) if (VEC_ELEM(*this,i)>max) {max=VEC_ELEM(*this,i); imax=i;} } /* Min index --------------------------------------------------------------- */ template void vT::min_index(int &imin) const { if (XSIZE(*this)==0) {imin=-1; return;} imin=0; T min=VEC_ELEM(*this,imin); FOR_ALL_ELEMENTS_IN_MATRIX1D(*this) if (VEC_ELEM(*this,i) void vT::compute_stats(double &avg, double &stddev, T &min_val, T &max_val, const matrix1D &corner1, const matrix1D &corner2) const { min_val=max_val=(*this)(corner1); matrix1D r(3); double N=0, sum=0, sum2=0; FOR_ALL_ELEMENTS_IN_MATRIX1D_BETWEEN(corner1,corner2) { sum+=(*this)(r); sum2+=(*this)(r)*(*this)(r); N++; if ((*this)(r)max_val) max_val=(*this)(r); } if (N!=0) { avg=sum/N; stddev=sqrt(sum2/N-avg*avg); } else {avg=stddev=0;} } /* Minimum and maximum in region ------------------------------------------- */ template void vT::compute_double_minmax(double &min_val, double &max_val, const matrix1D &corner1, const matrix1D &corner2) const { min_val=max_val=(*this)(corner1); matrix1D r(1); FOR_ALL_ELEMENTS_IN_MATRIX1D_BETWEEN(corner1,corner2) { if ((*this)(r)max_val) max_val=(*this)(r); } } /* Vector R2 and R3 -------------------------------------------------------- */ matrix1D vector_R2(double x, double y) { matrix1D result(2); VEC_ELEM(result,0)=x; VEC_ELEM(result,1)=y; return result; } matrix1D vector_R3(double x, double y, double z) { matrix1D result(3); VEC_ELEM(result,0)=x; VEC_ELEM(result,1)=y; VEC_ELEM(result,2)=z; return result; } matrix1D vector_R3(int x, int y, int z) { matrix1D result(3); VEC_ELEM(result,0)=x; VEC_ELEM(result,1)=y; VEC_ELEM(result,2)=z; return result; } /* Are orthogonal ---------------------------------------------------------- */ int are_orthogonal (matrix1D &v1, matrix1D &v2, matrix1D &v3) _THROW { if (XSIZE(v1)!=3 || XSIZE(v2)!=3 || XSIZE(v3)!=3) REPORT_ERROR(1002,"Are orthogonal: Some vector do not belong to R3"); try{ if (dot_product(v1,v2)!=0) return 0; if (dot_product(v2,v3)!=0) return 0; if (dot_product(v1,v3)!=0) return 0; } catch (Xmipp_error) { REPORT_ERROR(1007,"Are orthogonal: Vectors are not all of the same shape"); } return 1; } /* Are system? ------------------------------------------------------------- */ int are_system (matrix1D &v1, matrix1D &v2, matrix1D &v3) _THROW { matrix1D aux(3); if (XSIZE(v1)!=3 || XSIZE(v2)!=3 || XSIZE(v3)!=3) REPORT_ERROR(1002,"Are orthogonal: Some vector do not belong to R3"); aux=vector_product(v1,v2); if (aux!=v3) return 0; aux=vector_product(v2,v3); if (aux!=v1) return 0; aux=vector_product(v3,v1); if (aux!=v2) return 0; return 1; } /* Cut to common size ------------------------------------------------------ */ template void cut_to_common_size(vT &V1, vT &V2) { int x0=MAX(STARTINGX(V1) ,STARTINGX(V2)); int xF=MIN(FINISHINGX(V1),FINISHINGX(V2)); V1.window(x0,xF); V2.window(x0,xF); } /* Sort two vectors -------------------------------------------------------- */ template void sort_two_vectors(vT &v1, vT &v2) _THROW { T temp; if (XSIZE(v1)!=XSIZE(v2) || STARTINGX(v1)!=STARTINGX(v2)) REPORT_ERROR(1007, "sort_two_vectors: vectors are not of the same shape"); FOR_ALL_ELEMENTS_IN_MATRIX1D(v1) { temp=MIN(VEC_ELEM(v1,i),VEC_ELEM(v2,i)); VEC_ELEM(v2,i)=MAX(VEC_ELEM(v1,i),VEC_ELEM(v2,i)); VEC_ELEM(v1,i)=temp; } } /* Powell's optimizer ------------------------------------------------------ */ void Powell_optimizer(matrix1D &p, int i0, int n, double (*f)(double *x), double ftol, double &fret, int &iter, const matrix1D &steps, bool show) { double *xi=NULL; // Adapt indexes of p double *pptr=p.adapt_for_numerical_recipes(); double *auxpptr=pptr+(i0-1); // Form direction matrix ask_Tvector(xi,1,n*n); int ptr; for (int i=1,ptr=1; i<=n; i++) for (int j=1; j<=n; j++,ptr++) xi[ptr]=(i==j)?steps(i-1):0; // Optimize xi-=n; // This is because NR works with matrices // starting at [1,1] powell(auxpptr,xi,n,ftol,iter,fret,f, show); xi+=n; // Exit free_Tvector(xi,1,n*n); } /* Center of mass ---------------------------------------------------------- */ template void vT::center_of_mass(matrix1D ¢er, void * mask) { center.init_zeros(1); double mass=0; matrix1D *imask=(matrix1D *) mask; FOR_ALL_ELEMENTS_IN_MATRIX1D(*this) { if (imask==NULL || VEC_ELEM(*imask,i)) { XX(center)+=i*VEC_ELEM(*this,i); mass+=VEC_ELEM(*this,i); } } if (mass!=0) center/=mass; } /* The multidimensional arrays must be instantiated if we want them to compile into a library. All methods of all possible classes must be called. The only basic types considered are: short int double double complex For complex arrays not all methods are valid, so they are not in the library. */ /* ------------------------------------------------------------------------- */ /* INSTANTIATE */ /* ------------------------------------------------------------------------- */ template void instantiate_vector (matrix1D v) { matrix1D a; matrix1D r; matrix1D ir; double d; T Taux; // General functions for multidimensional arrays a==a; a!=a; a=1-a; a=a-1; a=a*a; a.print_stats(); a.compute_double_minmax(d,d); a.compute_double_minmax(d,d,r,r); a.compute_double_minmax(d,d,ir,ir); a.compute_stats(d,d,Taux,Taux,ir,ir); a.compute_max(); a.compute_min(); a.compute_avg(); a.compute_stddev(); a.range_adjust(0,1); a.statistics_adjust(0,1); a.effective_range(99); a.init_random(0,1); a.add_noise(0,1); a.threshold("abs_above",1,0); a.count_threshold("abs_above",1,0); a.center_of_mass(r); a.print_shape(); int size[3]; a.get_size(size); a.outside(r); a.isBorder(0); matrix1D pixel(1); a.isBorder(pixel); a.outside(0); a.intersects(a); a.intersects(r,r); a.isCorner(r); a.patch(a); cout << a; a.window(0,1); int imax; a.max_index(imax); a.min_index(imax); cut_to_common_size(a,a); // Specific for vectors matrix1D indx; a.init_linear(0,5); a.reverse(); a=a.sort(); indx=a.index_sort(); sort_two_vectors(a,a); } void instantiate_complex_vector () { matrix1D< complex > a; matrix1D r; // General functions for multidimensional arrays a=(complex)1.0-a; a=a-(complex)1.0; a=a*a; a.print_shape(); int size[3]; a.get_size(size); a.outside(r); a.outside(0); a.intersects(a); a.intersects(r,r); a.isCorner(r); a.patch(a); cout << a; a.window(0,1); cut_to_common_size(a,a); // Specific for vectors a.reverse(); } void instantiate1D() { matrix1D V0; instantiate_vector(V0); matrix1D V1; instantiate_vector(V1); matrix1D V2; instantiate_vector(V2); matrix1D V3; instantiate_vector(V3); matrix1D V4; instantiate_vector(V4); matrix1D V5; instantiate_vector(V5); } --=-vj1wM3x7pj21v3gE0Vq3--