/** @(#) matrixdouble/vectordouble.cpp */
#include "vectordouble.h"
// constructor, ~, friend-----------------------------------------------------
/** Construye un vector.
* Necesario un constructor sin parametros para crear arreglo de VectorDouble.
* @param theSize Taman~o del vector. Min: 1.
* @param theLogic 0: Logica 0. 1: Otro numero.
* @param cleanMe true: Limpia el vector con 0.0. false: no inicializa. */
VectorDouble :: VectorDouble(int theSize, int theLogic, bool cleanMe) {
//cout << "\n{ VectorDouble()";
v = NULL; // flags de creacion
resize(theSize, theLogic, cleanMe);
}
/** Constructor de copia. Copia dimension, logica y valores. */
VectorDouble :: VectorDouble(VectorDouble &source) {
int theSize = source.size;
int theLogic = source.logic;
//cout << "\n{ VectorDouble()";
v = NULL; // flags de creacion
resize(theSize, theLogic, false); // false se van a copiar
copyFrom(source);
}
/** Libera memoria. Para funciones fuera de clases, usar deleteArray().
* Para llamar al destructor theObject = NULL; */
VectorDouble :: ~VectorDouble() {
if (v != NULL) {
// Por algun error se tiene que atrapar
try {
//delete[] v;
} catch(exception e) {
cout << "~VectorDouble: " << e.what();
} catch(...) {
cout << "~VectorDouble: catch(...)";
}
// No usar v = NULL; // Asegurarse
}
}
/** Para ser usado por RectangularMatrix.
* Libera y asigna memoria. Se pierden los datos.
* Si es del mismo taman~o, no hace nada.
* Logica cero.*/
void VectorDouble :: resize(int theSize, int theLogic, bool cleanMe) {
/* algoritmo
resize(theSize); // 1.
logic = (theLogic ? 1 : 0); // 2. Necesita memoria
if (cleanMe) { // 3. Necesita logica
clean();
}
*/
int validSize = ((theSize > 0) ? theSize : 1);
// Primero se determina si se creo luego la validacion
// Si no esta creada, valida y reservar
if (v == NULL) {
v = new double[validSize];
// Si ya estaba creada y es del mismo taman~o, no hace nada
// Compara size anterior con el nuevo
} else if (size != validSize) {
// julio no lo hace
delete[] v; // 1o.
v = new double[validSize];
}
size = validSize;
logic = (theLogic ? 1 : 0); // 2. Necesita memoria
if (cleanMe) { // 3. Necesita logica
clean();
}
}
/** Asume dimensiones y logica iguales. */
// operator=
void VectorDouble :: copyFrom(VectorDouble &source) {
int i;
for (i = 0; i < size; ++i) {
v[i] = source.v[i];
}
}
// operator-------------------------------------------------------------------
/**
* Muestra el vector en la consola. No muestra las dimensiones.
* @param outStr Destino del flujo.
* @param vector Vector a mostrar.
* @return Vector
*/
ostream & operator<<(ostream &outStr, const VectorDouble &theVector) {
int i;
// Ajuste a la logica
int first = theVector.getLogic();
int last = theVector.getSize() + first - 1;
double value;
for (i = first; i <= last; ++i) {
value = theVector[i]; // error si << theVector[i]
outStr << value << "\t";
}
return outStr;
}
/**
* Entrada de la matriz en la consola.
* Permite lectura de archivo de la dimension y de los datos.
* No se lee la logica de la matriz.
* Uso:
* 3 1.1 2.2 3.3
* Logica variable.
* @param outStr Destino del flujo.
* @param matrix Matriz a mostrar.
* @return Matriz
*/
istream & operator>>(istream &inpStr, VectorDouble &theVector) {
//matrix.read(); // sig inpStr >>
int i, n;
inpStr >> n;
// Ajuste a la logica
int first = theVector.getLogic();
int last = n + first - 1;
// Reserva memoria
theVector.resize(n, first, true); // logic
for (i = first; i <= last; ++i) {
inpStr >> theVector[i];
// error si inpStr >> matrix(i, j)
}
return inpStr;
}
/**
* Permite asignacion (operator=).
* Logica 0, (0, 0) apunta al matematico A 1,1 esq sup izq
* @param row Logica cero: [0, row -1]. Logica uno: [1, row]
* @param column Logica cero: [0, col -1]. Logica uno: [1, col]
*/
double & VectorDouble :: operator()(int index)
throw (invalid_argument) {
if ( (index < logic) || (index >= (size + logic)) ) {
throw invalid_argument("operator(): indices fuera de rango.");
}
return v[index - logic];
}
/** @version const Sino error por rvalue. */
const double & VectorDouble :: operator()(int index) const
throw (invalid_argument) {
if ( (index < logic) || (index >= (size + logic)) ) {
throw invalid_argument("operator(): indices fuera de rango.");
}
return v[index - logic];
}
/** No checa dimensiones, ni logica. */
VectorDouble VectorDouble :: operator+(VectorDouble &right) {
VectorDouble res(size, logic);
int i;
for (i = 0; i < size; ++i) {
res.v[i] = v[i] + right.v[i];
}
return res;
}
/** No checa dimensiones, ni logica. */
VectorDouble VectorDouble :: operator-(VectorDouble &right) {
VectorDouble res(size, logic);
int i;
for (i = 0; i < size; ++i) {
res.v[i] = v[i] - right.v[i];
}
return res;
}
/**
* Multiplica un vector por un escalar.
* Logica cero.
* @param product Donde se almacenara el resultado
*/
void VectorDouble :: multiply(double escalar, VectorDouble &product)
throw (invalid_argument) {
if (size != product.size) {
throw invalid_argument("multiply: Dimensiones diferentes.");
}
int i;
for (i = 0; i < size; ++i) {
product.v[i] = v[i] * escalar;
}
}
/** Multiplica un vector por un escalar.
* Logica cero.
* @return Nuevo Vector multiplicado. */
// error de memoria
// dep
VectorDouble VectorDouble :: operator*(double right) {
int i;
VectorDouble res(size, logic); // mismo tam y log
for (i = 0; i < size; ++i) {
res.v[i] = right * v[i];
}
return res; // No se puede regresar solo product
}
/**
* Multiplica un vector por un escalar.
* Logica cero.
* @return Nuevo Vector multiplicado.
*/
/* error de memoria
VectorDouble VectorDouble :: operator*(double right) {
//double VectorDouble :: multiply(double escalar) {
int i;
VectorDouble product(size, logic); // mismo tam y log
for (i = 0; i < size; ++i) {
product[i + logic] = right * v[i];
}
return VectorDouble(product); // No se puede regresar solo product
}
*/
/** @version Usar constantes. */
/*VectorDouble VectorDouble :: operator*(const double right) {
return (*this) * const_cast<double>(right);
}
*/
/**
* Multiplica un vector por otro vector. Producto interno.
* Checa dimensiones.
* Logica cero.
* @param y Segundo operador.
*/
double VectorDouble :: operator*(VectorDouble &right) throw(invalid_argument) {
//double VectorDouble :: multiply(double escalar) {
if (size != right.getSize()) {
throw invalid_argument(
"operator*: Las dimensiones son diferentes.");
}
int i;
double sum = 0.0;
// proceso con la actual, pero manejar distintas logicas
// caso logic=1, rightLogic=0
for (i = 0; i < size; ++i) {
sum += v[i] * right.v[i];
}
return sum;
}
/** Redimensiona el vector, copia los datos de right y las caracteristicas. */
const VectorDouble &VectorDouble :: operator=(const VectorDouble &right) {
if (size != right.getSize()) {
resize(right.getSize(), right.getLogic(), false);
//throw invalid_argument(
// "operator=: Las dimensiones son diferentes.");
}
if (&right != this) { // Revisa autoasignacion
copyFrom( const_cast<VectorDouble &>(right) );
}
return *this; // Permite x = y = z;
}
// public---------------------------------------------------------------------
/**
* Llena el vector con el valor indicado.
* Logica cero.
* @param value Cada elemento valdra este numero.
*/
void VectorDouble :: fill(double value) {
int i;
for (i = 0; i < size; ++i) {
v[i] = value;
}
}
/**
* Solo lee los elementos sin la dimension desde consola.
* Logica cero.
*/
void VectorDouble :: readElements() {
int i;
for (i = 0; i < size; ++i) {
cin >> v[i];
}
}
/** @return true: Hay al menos un cero (menor que MACHEPS).
* false: todos los elementos distintos de cero.*/
bool VectorDouble :: hasAlmostZeros() {
int i;
for (i = 0; i < size; ++i) {
if ( fabs(v[i]) < MACHEPS) {
return true;
}
}
return false;
}
/** Usado por RectangularMatrix::deleteRectangular.
* Cuando no funciona el destructor (~). */
void VectorDouble :: deleteArray() {
if (v != NULL) {
delete[] v;
}
v = NULL;
}
/** @param digits Precision. Numero de digitos.
Logica cero. */
void VectorDouble :: print(int digits) {
int i;
cout << setiosflags(ios::fixed) << setprecision(digits);
for (i = 0; i < size; ++i) {
cout << v[i] << "\t";
}
cout << "\t";
}
/** @param digits Precision. Numero de digitos.
Logica cero. */
string VectorDouble :: toString(int digits, char separator) {
ostringstream outStr;
int i;
outStr << setiosflags(ios::fixed) << setprecision(digits);
for (i = 0; i < size; ++i) {
outStr << v[i] << separator;
}
return outStr.str();
}
/** Calcula la inversa. 1/v[i]. */
VectorDouble VectorDouble :: inverse() {
VectorDouble res(size, logic);
int i;
for (i = 0; i < size; ++i) {
res.v[i] = 1.0 / v[i];
}
return res;
}
/** Calcula el producto interno. Vector renglon * Vector columna = escalar.
* No checa dimensiones. */
// // ver operator *
/*double VectorDouble :: innerProduct(VectorDouble &right) {
double sum = 0.0;
int i;
for (i = 0; i < size; ++i) {
sum += v[i] * right[i];
}
return sum;
}
*/
/** Toma el vector actual como si fuera una matriz diagonal.
* [ v[0] 0 ][r[0]]
* [ v[2] ][r[1]]
* [ ... ][...]
* [ 0 v[n] ][r[n]]
* Se multiplica la matriz diagonal por un vector columna,
resultando un vector columna. */
VectorDouble VectorDouble :: diagonalProduct(VectorDouble &right) {
VectorDouble res(size, logic);
int i;
for (i = 0; i < size; ++i) {
res.v[i] = v[i] * right.v[i];
}
return res;
}
/** Calcula la norma superior o infinita.
* || v ||oo l(infinity) sup norm */
double VectorDouble :: infinityNorm() {
double mx = v[0];
int i;
// i = 1 se asigna el primero como maximo
for (i = 1; i < size; ++i) {
if (v[i] > mx) {
mx = v[i];
}
}
return mx;
}
/** Calcula la norma del minimo residuo absoluto
* || v ||1 l1 LAR norm (Least Absolute Residual) */
// No usar generalNorm por velocidad
double VectorDouble :: norm1() {
double sum = 0.0;
int i;
for (i = 0; i < size; ++i) {
sum += fabs(v[i]);
}
return sum;
}
/** Calcula la norma del minimo residuo absoluto
* || v ||2 l2 Euclidean or least-squares norm */
// No usar generalNorm por velocidad
double VectorDouble :: norm2() {
double sum = 0.0;
int i;
for (i = 0; i < size; ++i) {
sum += v[i] * v[i];
}
return sqrt(sum);
}
/** Calcula la norma del minimo residuo absoluto
* || v ||p lp General norm
* @param p p=[1, infinity) No lo checa. */
double VectorDouble :: generalNorm(double p) {
double sum = 0.0;
int i;
for (i = 0; i < size; ++i) {
sum += pow( fabs(v[i]), p);
}
return pow(sum, 1.0/p);
}
/** Lee de un archivo dos columnas. Ignora las cadenas. Ejemplo:
<pre>
Dim 3
Tiempo Poblacion
0.00 2.50000000000
0.25 4.46185192849
0.50 7.09001693225
</pre> */
//static
void VectorDouble :: read2ColumnsFromFile(string pathFile,
VectorDouble &v1, VectorDouble &v2)
throw (runtime_error) {
ifstream inFile(pathFile.data(), ios::in);
if ( !inFile ) { // !sobrecargado
throw runtime_error(
"graphDoglegFile: Error de apertura de archivo.");
}
string dummy;
int i, n;
// Dim
inFile >> dummy >> n;
v1.resize(n, 1, false);
v2.resize(n, 1, false);
// Titulo1 Titulo2
inFile >> dummy >> dummy;
for (i = 0; i < n; ++i) {
inFile >> v1.v[i] >> v2.v[i];
}
}
// Sin ajustar
/**
* compute norm1 for V[n] (vector)
* Logica uno.
*/
/*double norm1( double *V, int n) {
double sum = 0.0;
for (int i=1; i<=n; ++i)
sum += fabs(V[i]);
return sum;
}
*/
// Sin ajustar
/**
* Solve Rx = e for x .
* Logica uno.
*/
/*void solveRxe( double **R, double *x, int n) {
double *p = new double[n + 1];
double *pm = new double[n + 1];
double xp;
double xm;
double temp;
double tempm;
x[1] = 1.0/R[1][1];
for (int i=2; i<=n; ++i)
p[i] = x[1] * R[1][i];
for (int j=2; j<=n; ++j) {
xp = (1.0-p[j])/R[j][j];
xm = (-1.0-p[j])/R[j][j];
temp = fabs(xp);
tempm = fabs(xm);
for (int i=j+1; i<=n; ++i) {
pm[i] = p[i] + xm*R[j][i];
tempm += fabs(pm[i])/fabs(R[i][i]);
p[i] += xp *R[j][i];
temp += fabs(p[i])/fabs(R[i][i]);
}
if (temp >= tempm)
x[j] = xp;
else {
x[j] = xm;
for (int i=j+1; i<=n; ++i)
p[i] = pm[i];
}
}
delete []p;
delete []pm;
}
*/
// Sin ajustar
/**
* Solve Ry = x for y
* Logica uno.
*/
/*
void solveRyx(double **R, double *y, double *x, int n) {
double sum;
for (int i=n; i>=1; --i) {
sum = 0.0;
for (int j=i+1; j<=n; ++j)
sum += R[i][j] * y[j];
y[i] = (x[i] - sum)/R[i][i];
}
}
*/
// Sin ajustar
/**
* computes estimation for condition number k1(R)
* Logica uno.
*/
/*
double condest( double **R, int n ) {
double *x = new double[n + 1];
double *y = new double[n + 1];
double est;
est = norm1(R, n, n);
solveRxe(R, x, n);
solveRyx(R, y, x, n);
est *= norm1(y, n)/norm1(x, n);
delete []x;
delete []y;
return est;
}
*/
// private--------------------------------------------------------------------
// Fin------------------------------------------------------------------------