* Please, help with a weir error on templates.
@ 2004-03-17 15:34 Jose Roman Bilbao
[not found] ` <200403180007.51842.alexjdam@yahoo.com.br>
0 siblings, 1 reply; 2+ messages in thread
From: Jose Roman Bilbao @ 2004-03-17 15:34 UTC (permalink / raw)
To: gcc-help
Hi all,
I am not a member of this list, therefore I will be grateful if you CC
me a copy of your responses.
I have the next problem with a source file:
When compiler reaches the next code:
template <>
void core_array_by_scalar< complex<double> >(const maTC &op1,
const complex<double> &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
Well, I have included:
#include <complex>
Some time ago it was compiled with complex.h and it did works. Now it
does if I change it back. The most strage thing (for me) is that this
part of code compiles well with <complex> in the same file:
bool operator == (const ma< complex<double> > &op1, const ma<
complex<double> > &op2) {
On the other hand, if I change <complex.h> to something similat to that
found on complex.h so:
#include <complex>
using std::complex;
typedef complex<float> float_complex;
typedef complex<double> double_complex;
typedef complex<long double> long_double_complex;
It does compile perfectly... I am puzzled. Can you explain me which is
the error please?.
Thank you very much,
Roman
^ permalink raw reply [flat|nested] 2+ messages in thread
* Re: Please, help with a weir error on templates.
[not found] ` <200403180007.51842.alexjdam@yahoo.com.br>
@ 2004-03-18 15:19 ` Jose Roman Bilbao
0 siblings, 0 replies; 2+ messages in thread
From: Jose Roman Bilbao @ 2004-03-18 15:19 UTC (permalink / raw)
To: Alex J. Dam; +Cc: gcc-help
[-- Attachment #1: Type: text/plain, Size: 1168 bytes --]
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<double> >(const maTC &op1,
> > const complex<double> &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
[-- Attachment #2: MultidimBasic.inc --]
[-- Type: text/plain, Size: 16715 bytes --]
/***************************************************************************
*
* 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 <complex>
#define mi MULTIDIM_ELEM(*this,i)
#define msize MULTIDIM_SIZE(*this)
/* Print stats ------------------------------------------------------------- */
template <class T>
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 <class T>
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 <class T>
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<min_val) min_val=mi;
return min_val;
}
/* Compute min-max --------------------------------------------------------- */
template <class T>
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 (mi<min) min=(double)mi;
if (mi>max) max=(double)mi;
}
}
template <class T>
void maT::compute_double_minmax(double &min, double &max,
const matrix1D<int> &corner1, const matrix1D<int> &corner2) const {
matrix1D<double> dcorner1, dcorner2;
type_cast(corner1,dcorner1);
type_cast(corner2,dcorner2);
compute_double_minmax(min,max,dcorner1,dcorner2);
}
template <class T>
void maT::compute_stats(double &avg, double &stddev, T &min_val, T &max_val,
const matrix1D<int> &corner1, const matrix1D<int> &corner2) const {
matrix1D<double> dcorner1, dcorner2;
type_cast(corner1,dcorner1);
type_cast(corner2,dcorner2);
compute_stats(avg,stddev,min_val,max_val,dcorner1,dcorner2);
}
/* Compute avg ------------------------------------------------------------- */
template <class T>
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 <class T>
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 <class T>
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 (mi<min) min=mi;
}
avg /=msize;
if (msize>1) {
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 <class T>
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 <class T>
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 <class T>
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 <class T>
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 <class T>
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 <class T>
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=SGN(mi)*b; break;
case 3: if (mi>a) mi=b; break;
case 4: if (mi<a) mi=b; break;
case 5: if (mi<a) mi=a;
else if (mi>b) mi=b; break;
}
}
/* Count with threshold ---------------------------------------------------- */
template <class T>
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 3: if (mi>a) retval++; break;
case 4: if (mi<a) retval++; break;
case 5: if (mi>=a && mi<=b) retval++; break;
}
return retval;
}
/* Equality for normal data types ------------------------------------------ */
template <class T>
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<double> > &op1, const ma< complex<double> > &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 <class T>
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<double> >(const maTC &op1,
const complex<double> &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 <class T>
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<double> >(const complex<double> &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 <class T>
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<double> >(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;
}
}
[-- Attachment #3: xmippMatrices1D.cc --]
[-- Type: text/x-c++, Size: 17192 bytes --]
/***************************************************************************
*
* 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 <stdio.h>
/* ************************************************************************* */
/* IMPLEMENTATIONS */
/* ************************************************************************* */
#define maT matrix1D<T>
#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<T>, but
// we have "redirected" maT to matrix1D<T>, this way they will do OK in
// this vector library
/* Print shape ------------------------------------------------------------- */
template <class T>
void vT::print_shape(ostream &out) const {
out << "Size: " << XSIZE(*this)
<< "i=[" << STARTINGX(*this) << ".." << FINISHINGX(*this) << "]";
}
/* Get size--- ------------------------------------------------------------- */
template <class T>
void vT::get_size(int *size) const
{size[0]=xdim; size[1]=1; size[2]=1;}
/* Outside ----------------------------------------------------------------- */
template <class T>
bool vT::outside(const matrix1D<double> &v) const _THROW {
if (XSIZE(v)<1)
REPORT_ERROR(1,"Outside: index vector has got not enough components");
return (XX(v)<STARTINGX(*this) || XX(v)>FINISHINGX(*this));
}
template <class T>
bool vT::outside(int i) const {
return (i<STARTINGX(*this) || i>FINISHINGX(*this));
}
/* Intersects -------------------------------------------------------------- */
template <class T>
bool vT::intersects(const vT &m) const
{return intersects(STARTINGX(m), XSIZE(m)-1);}
template <class T>
bool vT::intersects(const matrix1D<double> &corner1,
const matrix1D<double> &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 <class T>
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 <class T>
bool vT::isCorner(const matrix1D<double> &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 <class T>
bool vT::isBorder(const matrix1D<int> &v) _THROW
{
if (XSIZE(v)<1)
REPORT_ERROR(1,"isBorder: index vector has got not enough components");
return isBorder(XX(v));
}
template <class T>
bool vT::isBorder(int i)
{
return (i==STARTINGX(*this) || i==FINISHINGX(*this));
}
/* Patch ------------------------------------------------------------------- */
template <class T>
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 <class T>
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<double> >& 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 <class T>
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<steps; i++)
VEC_ELEM(*this,i)=(T) ((double)minF+slope*i);
}
}
/* Reverse ----------------------------------------------------------------- */
template <class T>
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 <class T>
vT vT::sort() const {
vT temp;
matrix1D<double> 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 <class T>
matrix1D<int> vT::index_sort() const {
matrix1D<int> indx;
matrix1D<double> 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 <class T>
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 <class T>
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 <class T>
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)<min) {min=VEC_ELEM(*this,i); imin=i;}
}
/* Statistics in region ---------------------------------------------------- */
template <class T>
void vT::compute_stats(double &avg, double &stddev, T &min_val, T &max_val,
const matrix1D<double> &corner1, const matrix1D<double> &corner2) const {
min_val=max_val=(*this)(corner1);
matrix1D<double> 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)<min_val) min_val=(*this)(r);
else 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 <class T>
void vT::compute_double_minmax(double &min_val, double &max_val,
const matrix1D<double> &corner1, const matrix1D<double> &corner2) const {
min_val=max_val=(*this)(corner1);
matrix1D<double> r(1);
FOR_ALL_ELEMENTS_IN_MATRIX1D_BETWEEN(corner1,corner2) {
if ((*this)(r)<min_val) min_val=(*this)(r);
else if ((*this)(r)>max_val) max_val=(*this)(r);
}
}
/* Vector R2 and R3 -------------------------------------------------------- */
matrix1D<double> vector_R2(double x, double y) {
matrix1D<double> result(2);
VEC_ELEM(result,0)=x;
VEC_ELEM(result,1)=y;
return result;
}
matrix1D<double> vector_R3(double x, double y, double z) {
matrix1D<double> result(3);
VEC_ELEM(result,0)=x;
VEC_ELEM(result,1)=y;
VEC_ELEM(result,2)=z;
return result;
}
matrix1D<int> vector_R3(int x, int y, int z) {
matrix1D<int> 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<double> &v1, matrix1D<double> &v2,
matrix1D<double> &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<double> &v1, matrix1D<double> &v2,
matrix1D<double> &v3) _THROW {
matrix1D<double> 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 <class T>
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 <class T>
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<double> &p, int i0, int n,
double (*f)(double *x), double ftol, double &fret,
int &iter, const matrix1D<double> &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 <class T>
void vT::center_of_mass(matrix1D<double> ¢er, void * mask) {
center.init_zeros(1);
double mass=0;
matrix1D<int> *imask=(matrix1D<int> *) 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<double>
For complex<double> arrays not all methods are valid, so they are not in
the library.
*/
/* ------------------------------------------------------------------------- */
/* INSTANTIATE */
/* ------------------------------------------------------------------------- */
template <class T>
void instantiate_vector (matrix1D<T> v) {
matrix1D<T> a;
matrix1D<double> r;
matrix1D<int> 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<int> 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<int> 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<double> > a;
matrix1D<double> r;
// General functions for multidimensional arrays
a=(complex<double>)1.0-a;
a=a-(complex<double>)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<char> V0; instantiate_vector(V0);
matrix1D<short> V1; instantiate_vector(V1);
matrix1D<int> V2; instantiate_vector(V2);
matrix1D<float> V3; instantiate_vector(V3);
matrix1D<double> V4; instantiate_vector(V4);
matrix1D<long double> V5; instantiate_vector(V5);
}
[-- Attachment #4: xmippMatrices1D.cc --]
[-- Type: text/x-c++, Size: 17192 bytes --]
/***************************************************************************
*
* 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 <stdio.h>
/* ************************************************************************* */
/* IMPLEMENTATIONS */
/* ************************************************************************* */
#define maT matrix1D<T>
#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<T>, but
// we have "redirected" maT to matrix1D<T>, this way they will do OK in
// this vector library
/* Print shape ------------------------------------------------------------- */
template <class T>
void vT::print_shape(ostream &out) const {
out << "Size: " << XSIZE(*this)
<< "i=[" << STARTINGX(*this) << ".." << FINISHINGX(*this) << "]";
}
/* Get size--- ------------------------------------------------------------- */
template <class T>
void vT::get_size(int *size) const
{size[0]=xdim; size[1]=1; size[2]=1;}
/* Outside ----------------------------------------------------------------- */
template <class T>
bool vT::outside(const matrix1D<double> &v) const _THROW {
if (XSIZE(v)<1)
REPORT_ERROR(1,"Outside: index vector has got not enough components");
return (XX(v)<STARTINGX(*this) || XX(v)>FINISHINGX(*this));
}
template <class T>
bool vT::outside(int i) const {
return (i<STARTINGX(*this) || i>FINISHINGX(*this));
}
/* Intersects -------------------------------------------------------------- */
template <class T>
bool vT::intersects(const vT &m) const
{return intersects(STARTINGX(m), XSIZE(m)-1);}
template <class T>
bool vT::intersects(const matrix1D<double> &corner1,
const matrix1D<double> &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 <class T>
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 <class T>
bool vT::isCorner(const matrix1D<double> &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 <class T>
bool vT::isBorder(const matrix1D<int> &v) _THROW
{
if (XSIZE(v)<1)
REPORT_ERROR(1,"isBorder: index vector has got not enough components");
return isBorder(XX(v));
}
template <class T>
bool vT::isBorder(int i)
{
return (i==STARTINGX(*this) || i==FINISHINGX(*this));
}
/* Patch ------------------------------------------------------------------- */
template <class T>
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 <class T>
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<double> >& 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 <class T>
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<steps; i++)
VEC_ELEM(*this,i)=(T) ((double)minF+slope*i);
}
}
/* Reverse ----------------------------------------------------------------- */
template <class T>
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 <class T>
vT vT::sort() const {
vT temp;
matrix1D<double> 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 <class T>
matrix1D<int> vT::index_sort() const {
matrix1D<int> indx;
matrix1D<double> 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 <class T>
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 <class T>
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 <class T>
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)<min) {min=VEC_ELEM(*this,i); imin=i;}
}
/* Statistics in region ---------------------------------------------------- */
template <class T>
void vT::compute_stats(double &avg, double &stddev, T &min_val, T &max_val,
const matrix1D<double> &corner1, const matrix1D<double> &corner2) const {
min_val=max_val=(*this)(corner1);
matrix1D<double> 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)<min_val) min_val=(*this)(r);
else 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 <class T>
void vT::compute_double_minmax(double &min_val, double &max_val,
const matrix1D<double> &corner1, const matrix1D<double> &corner2) const {
min_val=max_val=(*this)(corner1);
matrix1D<double> r(1);
FOR_ALL_ELEMENTS_IN_MATRIX1D_BETWEEN(corner1,corner2) {
if ((*this)(r)<min_val) min_val=(*this)(r);
else if ((*this)(r)>max_val) max_val=(*this)(r);
}
}
/* Vector R2 and R3 -------------------------------------------------------- */
matrix1D<double> vector_R2(double x, double y) {
matrix1D<double> result(2);
VEC_ELEM(result,0)=x;
VEC_ELEM(result,1)=y;
return result;
}
matrix1D<double> vector_R3(double x, double y, double z) {
matrix1D<double> result(3);
VEC_ELEM(result,0)=x;
VEC_ELEM(result,1)=y;
VEC_ELEM(result,2)=z;
return result;
}
matrix1D<int> vector_R3(int x, int y, int z) {
matrix1D<int> 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<double> &v1, matrix1D<double> &v2,
matrix1D<double> &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<double> &v1, matrix1D<double> &v2,
matrix1D<double> &v3) _THROW {
matrix1D<double> 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 <class T>
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 <class T>
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<double> &p, int i0, int n,
double (*f)(double *x), double ftol, double &fret,
int &iter, const matrix1D<double> &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 <class T>
void vT::center_of_mass(matrix1D<double> ¢er, void * mask) {
center.init_zeros(1);
double mass=0;
matrix1D<int> *imask=(matrix1D<int> *) 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<double>
For complex<double> arrays not all methods are valid, so they are not in
the library.
*/
/* ------------------------------------------------------------------------- */
/* INSTANTIATE */
/* ------------------------------------------------------------------------- */
template <class T>
void instantiate_vector (matrix1D<T> v) {
matrix1D<T> a;
matrix1D<double> r;
matrix1D<int> 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<int> 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<int> 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<double> > a;
matrix1D<double> r;
// General functions for multidimensional arrays
a=(complex<double>)1.0-a;
a=a-(complex<double>)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<char> V0; instantiate_vector(V0);
matrix1D<short> V1; instantiate_vector(V1);
matrix1D<int> V2; instantiate_vector(V2);
matrix1D<float> V3; instantiate_vector(V3);
matrix1D<double> V4; instantiate_vector(V4);
matrix1D<long double> V5; instantiate_vector(V5);
}
^ permalink raw reply [flat|nested] 2+ messages in thread
end of thread, other threads:[~2004-03-18 9:15 UTC | newest]
Thread overview: 2+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2004-03-17 15:34 Please, help with a weir error on templates Jose Roman Bilbao
[not found] ` <200403180007.51842.alexjdam@yahoo.com.br>
2004-03-18 15:19 ` Jose Roman Bilbao
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).