diff --git a/expui/#BasisFactory.cc# b/expui/#BasisFactory.cc# new file mode 100644 index 000000000..8c5d01692 --- /dev/null +++ b/expui/#BasisFactory.cc# @@ -0,0 +1,380 @@ +p -#include + +#include +#include +#include +#include +#include +#include + +#ifdef HAVE_FE_ENABLE +#include +#endif + +namespace BasisClasses +{ + std::map Basis::coordLabels = + { {Basis::Coord::Spherical, "Spherical" }, + {Basis::Coord::Cylindrical, "Cylindrical"}, + {Basis::Coord::Cartesian, "Cartesian" }, + {Basis::Coord::None, "None" } }; + + Basis::Basis(const YAML::Node& CONF, const std::string& id) + { + // Class name + // + name = id; + + // Copy the YAML config + // + node = CONF; + + // Complete the initialization + // + initialize(); + } + + Basis::Basis(const std::string& confstr, const std::string& id) + { + // Assign class name + // + name = id; + + // Read the YAML from a string + // + try { + node = YAML::Load(confstr); + } + catch (const std::runtime_error& error) { + std::cout << "Basis constructor: found a problem in the YAML config" + << std::endl; + throw; + } + + // Complete the initialization + // + initialize(); + } + + void Basis::initialize() + { +#ifdef HAVE_FE_ENABLE + // Flag invalid FP results only, such as 0/0 or infinity - infinity + // or sqrt(-1). + // + // feenableexcept(FE_INVALID); +#endif + + // Check whether MPI is initialized + // + int flag; + MPI_Initialized(&flag); + if (flag) use_mpi = true; + else use_mpi = false; + + // Fall back sanity (works for me but this needs to be fixed + // generally) + // + if (use_mpi) { + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + } + + // Parameters for force + // + try { + conf = node["parameters"]; + } + catch (YAML::Exception & error) { + if (myid==0) std::cout << "Error parsing Basis parameters for <" + << name << ">: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << node << std::endl + << std::string(60, '-') << std::endl; + + throw std::runtime_error("Basis: error parsing YAML"); + } + + // Set coefficient center to zero by default + // + coefctr = {0.0, 0.0, 0.0}; + + // Default null coordinate type + // + coordinates = Coord::None; + } + + Basis::Coord Basis::parseFieldType(std::string coord_type) + { + // Find coordinate type + // + std::transform(coord_type.begin(), coord_type.end(), coord_type.begin(), + [](unsigned char c){ return std::tolower(c); }); + + Coord ctype; + + if (coord_type.find("cyl") != std::string::npos) { + ctype = Coord::Cylindrical; + } else if (coord_type.find("cart") != std::string::npos) { + ctype = Coord::Cartesian; + } else if (coord_type.find("none") != std::string::npos) { + ctype = Coord::None; + } else { + ctype = Coord::Spherical; + } + + return ctype; + } + + std::shared_ptr Basis::factory_string(const std::string& conf) + { + YAML::Node node; + + try { + // Read the YAML from a string + // + node = YAML::Load(conf); + } + catch (const std::runtime_error& error) { + std::cout << "Basis constructor: found a problem in the YAML config" + << std::endl; + throw; + } + + // Complete the initialization + // + return factory_initialize(node); + } + + std::shared_ptr Basis::factory(const YAML::Node& conf) + { + return factory_initialize(conf); + } + + + std::shared_ptr Basis::factory_initialize(const YAML::Node& conf) + { + std::shared_ptr basis; + std::string name; + + // Load parameters from YAML configuration node + try { + name = conf["id"].as(); + } + catch (YAML::Exception & error) { + if (myid==0) std::cout << "Error parsing force id in BasisFactory" + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + + throw std::runtime_error("BasisFactory: error parsing YAML"); + } + + try { + if ( !name.compare("sphereSL") ) { + basis = std::make_shared(conf); + } + else if ( !name.compare("bessel") ) { + basis = std::make_shared(conf); + } + else if ( !name.compare("cylinder") ) { + basis = std::make_shared(conf); + } + else if ( !name.compare("flatdisk") ) { + basis = std::make_shared(conf); + } + else if ( !name.compare("CBDisk") ) { + basis = std::make_shared(conf); + } + else if ( !name.compare("slabSL") ) { + basis = std::make_shared(conf); + } + else if ( !name.compare("cube") ) { + basis = std::make_shared(conf); + } + else if ( !name.compare("field") ) { + basis = std::make_shared(conf); + } + else if ( !name.compare("velocity") ) { + basis = std::make_shared(conf); + } + else { + std::string msg("I don't know about the basis named: "); + msg += name; + msg += ". Known types are currently 'sphereSL', 'cylinder', 'flatdisk', 'slabSL', 'field', and 'velocity'"; + throw std::runtime_error(msg); + } + } + catch (std::exception& e) { + std::cout << "Error in BasisFactory constructor: " << e.what() << std::endl; + throw; // Rethrow the exception? + } + + return basis; + } + + std::vector + Basis::operator()(double x1, double x2, double x3, const Coord ctype) + { + if (ctype==Coord::Spherical) + return sph_eval(x1, x2, x3); + else if (ctype==Coord::Cylindrical) + return cyl_eval(x1, x2, x3); + else if (ctype==Coord::Cartesian) + return crt_eval(x1, x2, x3); + else { + return sph_eval(x1, x2, x3); + }; + } + + std::vector Basis::getFields(double x, double y, double z) + { + return crt_eval(x, y, z); + } + + std::tuple, Eigen::VectorXd> + Basis::getFieldsCoefs + (double x, double y, double z, std::shared_ptr coefs) + { + // Python dictonary for return + std::map ret; + + // Times for the coefficients + auto times = coefs->Times(); + + // Initialize the dictionary/map + auto fields = getFieldLabels(coordinates); + for (auto s : fields) ret[s].resize(times.size()); + + // Make the return dictionary of arrays + for (int i=0; igetCoefStruct(times[i])); + // The field evaluation + auto v = crt_eval(x, y, z); + // Pack the fields into the dictionary + for (int j=0; j(times.data(), times.size()); + + // Return the dictionary and the time array + return {ret, T}; + } + + // Generate coefficients from the accumulated array values + CoefClasses::CoefStrPtr Basis::makeFromArray(double time) + { + make_coefs(); + load_coefs(coefret, time); + return coefret; + } + + // Generate coefficients from a phase-space table + // + CoefClasses::CoefStrPtr Basis::createFromArray + (Eigen::VectorXd& m, RowMatrixXd& p, double time, std::vector ctr, + bool roundrobin, bool posvelrows) + { + initFromArray(ctr); + addFromArray(m, p, roundrobin, posvelrows); + return makeFromArray(time); + } + + void Basis::setNonInertial(int N, Eigen::VectorXd& t, Eigen::MatrixXd& pos) + { + Naccel = N; + t_accel = t; + p_accel = pos; + } + + void Basis::setNonInertial(int N, const std::string& orient) + { + std::ifstream in(orient); + + if (not in) { + throw std::runtime_error("Cannot open Orient file with centering data: " + orient); + } + + const int cbufsiz = 16384; + std::unique_ptr cbuf(new char [cbufsiz]); + + // Look for data and write it while + // accumlating data for averaging + Eigen::Vector3d testread; + double time, dummy; + + std::vector times; + std::vector centers; + + while (in) { + + in.getline(cbuf.get(), cbufsiz); + if (in.rdstate() & (ios::failbit | ios::eofbit)) break; + + // Skip comment lines + // + if (cbuf[0] == '#') continue; + + std::istringstream line(cbuf.get()); + + // Read until current time is reached + line >> time; // + line >> dummy; + line >> dummy; + + bool allRead = true; + for (int i=0; i<8; i++) { + if (line.eof()) allRead = false; + for (int k; k<3; k++) line >> testread(k); + } + if (allRead) { + times.push_back(time); + centers.push_back(testread); + } + } + + // Repack data + Naccel = N; + t_accel.resize(times.size()); + p_accel.resize(times.size(), 3); + for (int i=0; i t_accel(t_accel.size()-1)) { + std::cout << "ERROR: " << time << " is outside of range of the non-inertial DB" + << std::endl; + ret.setZero(); + return ret; + } + else { + int imin = 0; + int imax = std::lower_bound(t_accel.data(), t_accel.data()+t_accel.size(), time) - t_accel.data(); + if (imax > Naccel) imin = imax - Naccel; + + int num = imax - imin + 1; + Eigen::VectorXd t(num); + Eigen::MatrixXd p(num, 3); + for (int i=imin; i<=imax; i++) { + t(i-imin) = t_accel(i); + for (int k=0; k<3; k++) p(i-imin, k) = p_accel(i, k); + } + + for (int k=0; k<3; k++) + ret(k) = 2.0*std::get<0>(QuadLS(t, p.col(k)).coefs()); + } + return ret; + } + +} +// END namespace BasisClasses diff --git a/expui/#SvdSignChoice.cc# b/expui/#SvdSignChoice.cc# new file mode 100644 index 000000000..187df58e2 --- /dev/null +++ b/expui/#SvdSignChoice.cc# @@ -0,0 +1,287 @@ +#include + +namespace MSSA { + + // SVD with corrected signs + // + // Input + // ----- + // X is data matrix (constant) + // U is the matrix of left singular vectors + // S is the array of singular values (constant) + // V is the matrix of right singular vectors + // + // Output + // ------ + // U, V returned corrected, disambiguated signs + // + // Reference + // --------- + // Bro, R., Acar, E., & Kolda, T. G. (2008). Resolving the sign + // ambiguity in the singular value decomposition. Journal of + // Chemometrics: A Journal of the Chemometrics Society, 22(2), + // 135-140. + // + // URL: + // https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2007/076422.pdf + // + void SvdSignChoice + (const Eigen::MatrixXd& X, + Eigen::MatrixXd& U, const Eigen::VectorXd& S, Eigen::MatrixXd& V, + bool stride) + { + // SDV dimensions + int I = U.rows(); + int J = V.rows(); + int K = S.size(); + + // Dimensions + // ---------- + // X is a [I x J] data matrix + // U is a [I x K] matrix + // S is a [K x 1] vector (diagonal matrix) + // V is a [J x K] matrix + + // Sanity checks + if (U.cols() != K) + throw std::invalid_argument("SvdSignChoice: U has wrong dimensions"); + + if (V.cols() != K) + throw std::invalid_argument("SvdSignChoice: V has wrong dimensions"); + + if (X.rows() != I || X.cols() != J) + throw std::invalid_argument("SvdSignChoice: X dimensions do not match SVD input"); + + // Check for stride + // + Eigen::MatrixXd colX, rowX; + if (stride) { + int minDim = std::min({I, J, K}); + int newDim = 2*minDim; + int Cstd = X.cols()/newDim, Rstd = X.rows()/newDim; + if (Cstd < 2 || Rstd < 2) + stride = false; + else { + colX.resize(I, newDim); + for (int j=0; j int + { + return (0.0 < val) - (val < 0.0); + }; + + // Determine and apply the sign correction + // +#pragma omp parallel for + for (int k=0; k int + { + return (0.0 < val) - (val < 0.0); + }; + + // Determine and apply the sign correction + // +#pragma omp parallel for + for (int k=0; k int + { + return (0.0 < val) - (val < 0.0); + }; + + // Determine and apply the sign correction + // +#pragma omp parallel for + for (int k=0; k getFieldLabels(void) { return getFieldLabels(coordinates); } + //! Get the basis expansion center + std::vector getCenter() { return coefctr; } }; using BasisPtr = std::shared_ptr; diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index f9237bc65..e6eab3a2d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3685,11 +3685,17 @@ namespace BasisClasses // auto basis = std::get<0>(mod); + // Get expansion center + // + auto ctr = basis->getCenter(); + // Get fields // int rows = accel.rows(); for (int n=0; ngetFields(ps(n, 0), ps(n, 1), ps(n, 2)); + auto v = basis->getFields(ps(n, 0) - ctr[0], + ps(n, 1) - ctr[1], + ps(n, 2) - ctr[2]); // First 6 fields are density and potential, follewed by acceleration for (int k=0; k<3; k++) accel(n, k) += v[6+k] - basis->pseudo(k); } @@ -3899,10 +3905,14 @@ namespace BasisClasses // int numT = floor( (tfinal - tinit)/h ); - // Compute output step + // Number of output steps. Will attempt to find the best stride... + // + int stride = std::ceil(static_cast(numT)/static_cast(nout)); + if (stride>1) numT = nout * stride; + + // Compute the interval-matching step // - nout = std::min(numT, nout); - double H = (tfinal - tinit)/nout; + h = (tfinal - tinit)/(numT - 1); // Return data // @@ -3934,9 +3944,11 @@ namespace BasisClasses for (int k=0; k<6; k++) ret(n, k, 0) = ps(n, k); double tnow = tinit; - for (int s=1, cnt=1; s= H*cnt-h*1.0e-8) { + if (s++ % stride == 0) { times(cnt) = tnow; for (int n=0; n +#include +#include + +#include +#include +#include + +// Verbose output for checking stack usage +#ifdef DEBUG +#include +#include +#endif + +namespace Field +{ + + FieldGenerator::FieldGenerator(const std::vector &time, + const std::vector &pmin, + const std::vector &pmax, + const std::vector &grid + ) : + times(time), pmin(pmin), pmax(pmax), grid(grid) + { + // Check whether MPI is initialized + // + int flag; + MPI_Initialized(&flag); + if (flag) use_mpi = true; + else use_mpi = false; + + + // Fall back sanity (works for me but this needs to be fixed + // generally) + // + if (use_mpi) { + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + } + } + + void FieldGenerator::check_times(CoefClasses::CoefsPtr coefs) + { + std::vector ctimes = coefs->Times(); + std::sort(ctimes.begin(), ctimes.end()); + + for (auto t : times) { + if (std::find(ctimes.begin(), ctimes.end(), t) == ctimes.end()) { + std::ostringstream sout; + sout << "FieldGenerator: requested time <" << t << "> " + << "not in DB" << std::endl; + throw std::runtime_error(sout.str()); + } + } + } + + std::map> + FieldGenerator::lines + (BasisClasses::BasisPtr basis, CoefClasses::CoefsPtr coefs, + std::vector beg, std::vector end, int num) + { + // Check + // + if (beg.size()!=3 or end.size()!=3) { + throw std::runtime_error("FieldGenerator::lines: vectors beg and end must have rank 3"); + } + + if (num<1) { + throw std::runtime_error("FieldGenerator::lines: number of evaluation points must be > 0"); + } + + check_times(coefs); + + std::map> ret; + + std::vector labels = + {"x", "y", "z", "arc", + "potl", "potl m=0", + "rad force", "mer force", "azi force", + "dens", "dens m=0"}; + + + std::vector dd(3); + for (int k=0; k<3; k++) dd[k] = (end[k] - beg[k])/num; + double dlen = sqrt(dd[0]*dd[0] + dd[1]*dd[1] + dd[2]*dd[2]); + + for (auto T : times) { + + if (not coefs->getCoefStruct(T)) { + std::cout << "Could not find time=" << T << ", continuing" << std::endl; + continue; + } + + basis->set_coefs(coefs->getCoefStruct(T)); + + std::map frame; + for (auto label : labels) { + frame[label] = Eigen::VectorXf::Zero(num); + } + + double r, phi, costh; + double p0, p1, d0, d1, fr, ft, fp; + + int ncnt = 0; // Process counter for MPI + + for (int ncnt=0; ncntall_eval(r, costh, phi, d0, d1, p0, p1, fr, ft, fp); + + frame["x" ](ncnt) = x; + frame["y" ](ncnt) = y; + frame["z" ](ncnt) = z; + frame["arc" ](ncnt) = dlen*ncnt; + frame["potl" ](ncnt) = p0 + p1; + frame["potl m=0" ](ncnt) = p0; + frame["rad force"](ncnt) = fr; + frame["mer force"](ncnt) = ft; + frame["axi force"](ncnt) = fp; + frame["dens" ](ncnt) = d0 + d1; + frame["dens m=0" ](ncnt) = d0; + } + } + + if (use_mpi) { + for (auto & f : frame) { + if (myid==0) + MPI_Reduce(MPI_IN_PLACE, f.second.data(), f.second.size(), + MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD); + else + MPI_Reduce(f.second.data(), NULL, f.second.size(), + MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD); + } + } + + ret[T] = frame; + } + + return ret; + } + + void FieldGenerator::file_lines + (BasisClasses::BasisPtr basis, CoefClasses::CoefsPtr coefs, + std::vector beg, std::vector end, int num, + const std::string prefix, const std::string outdir) + { + auto db = lines(basis, coefs, beg, end, num); + + if (myid==0) { + + int icnt = 0; + + for (auto & frame : db) { + + // Create the file name and open + // + std::ostringstream sout; + sout << outdir << "/" << prefix << "_probe_" << icnt << ".txt"; + + std::ofstream out(sout.str()); + + if (out) { + + // Write file header + // + out << "# T=" << frame.first << std::endl; + + int cnt = 0; + for (auto u : frame.second) { + std::string label = u.first + " "; + if (cnt++==0) out << "#" << std::setw(15) << std::right << label; + else out << std::setw(16) << std::right << label; + } + out << std::endl; + for (int n=0; n"); + } + + icnt++; + } + } + } + + + std::map> + FieldGenerator::slices(BasisClasses::BasisPtr basis, + CoefClasses::CoefsPtr coefs) + { + // Check + // + check_times(coefs); + + std::map> ret; + + std::vector labels = + {"potl", "potl m=0", "rad force", "mer force", "azi force", + "dens", "dens m=0"}; + + // Find the first two non-zero indices + // + int i1=-1, i2=-1, i3=-1; + std::vector pos, del; + for (size_t i=0; i0) { + if (i1<0) { + i1 = i; + del.push_back((pmax[i1] - pmin[i1])/std::max(grid[i1]-1, 1)); + } + else if (i2<0) { + i2 = i; + del.push_back((pmax[i2] - pmin[i2])/std::max(grid[i2]-1, 1)); + } + } else { + del.push_back(0.0); + i3 = i; + } + } + + if (i1<0 or i2<0 or i3<0) + throw std::runtime_error("FieldGenerator::slices: bad grid specification"); + + int ncnt = 0; // Process counter for MPI + + std::map frame; + for (auto label : labels) { + frame[label].resize(grid[i1], grid[i2]); + } + + for (auto T : times) { + + if (ncnt++ % numprocs > 0) continue; + + if (not coefs->getCoefStruct(T)) { + std::cout << "Could not find time=" << T << ", continuing" << std::endl; + continue; + } + + basis->set_coefs(coefs->getCoefStruct(T)); + + double r, phi, costh; + double p0, p1, d0, d1, fr, ft, fp; + + int totpix = grid[i1] * grid[i2]; + +#pragma omp parallel for + for (int k=0; kall_eval(r, costh, phi, + d0, d1, p0, p1, fr, ft, fp); + + // Pack the frame structure + // + frame["potl" ](i, j) = p0 + p1; + frame["potl m=0" ](i, j) = p0; + frame["rad force"](i, j) = fr; + frame["mer force"](i, j) = ft; + frame["azi force"](i, j) = fp; + frame["dens" ](i, j) = d0 + d1; + frame["dens m=0" ](i, j) = d0; + } + + ret[T] = frame; + } + + if (use_mpi) { + + std::vector bf(9); + + for (int n=1; n in field map" + << std::endl; + } + + // Get the data + // + MPI_Recv(frame[s].data(), frame[s].size(), MPI_FLOAT, n, 106, + MPI_COMM_WORLD, &status); + } + + ret[T] = frame; + } + } + } + } + + return ret; + } + + void FieldGenerator::file_slices(BasisClasses::BasisPtr basis, + CoefClasses::CoefsPtr coefs, + const std::string prefix, + const std::string outdir) + { + auto db = slices(basis, coefs); + + if (myid==0) { + + // Find the first two non-zero indices + int i1=-1, i2=-1, i3=-1; + for (size_t i=0; i0) { + if (i1<0) i1 = i; + else if (i2<0) i2 = i; + } else i3 = i; + } + + int icnt = 0; + + for (auto & frame : db) { + + DataGrid datagrid(grid[i1], grid[i2], 1, + pmin[i1], pmax[i1], pmin[i2], pmax[i2], 0, 0); + + std::vector tmp(grid[i1]*grid[i2]); + + for (auto & v : frame.second) { + + for (int i=0; i>> + FieldGenerator::volumes(BasisClasses::BasisPtr basis, + CoefClasses::CoefsPtr coefs) + { + std::map>> ret; + + std::vector labels = + {"x", "y", "z", "arc", + "potl", "potl m=0", + "rad force", "mer force", "azi force", + "dens", "dens m=0"}; + + // Allocate frame storge + // + std::map> frame; + for (auto label : labels) { + frame[label].resize({grid[0]+1, grid[1]+1, grid[2]+1}); + } + + std::vector del = + { (pmax[0] - pmin[0])/std::max(grid[0]-1, 1), + (pmax[1] - pmin[1])/std::max(grid[1]-1, 1), + (pmax[2] - pmin[2])/std::max(grid[2]-1, 1) }; + + int ncnt = 0; // Process counter for MPI + + for (auto T : times) { + + if (ncnt++ % numprocs > 0) continue; + + basis->set_coefs(coefs->getCoefStruct(T)); + + double r, phi, costh; + double p0, p1, d0, d1, fr, ft, fp; + + int totpix = grid[0] * grid[1] * grid[2]; + +#pragma omp parallel for + for (int n=0; nall_eval(r, costh, phi, d0, d1, p0, p1, fr, ft, fp); + + // Pack the frame structure + // + frame["potl" ](i, j, k) = p0 + p1; + frame["potl m=0" ](i, j, k) = p0; + frame["rad force"](i, j, k) = fr; + frame["mer force"](i, j, k) = ft; + frame["azi force"](i, j, k) = fp; + frame["dens" ](i, j, k) = d0 + d1; + frame["dens m=0" ](i, j, k) = d0; + } + + ret[T] = frame; + +#ifdef DEBUG + if (myid==0) { + rusage usage; + int err = getrusage(RUSAGE_SELF, &usage); + std::cout << "volumes: T=" << std::setw(8) << std::fixed<< T + << " Size=" << std::setw(8) << usage.ru_maxrss/1024/1024 + << std::endl; + } +#endif + + } + + if (use_mpi) { + + std::vector bf(9); + + for (int n=1; n in field map" + << std::endl; + } + + // Get the data + // + MPI_Recv(frame[s].data(), frame[s].size(), MPI_FLOAT, n, 106, + MPI_COMM_WORLD, &status); + } + + ret[T] = frame; +#ifdef DEBUG + rusage usage; + int err = getrusage(RUSAGE_SELF, &usage); + std::cout << "volumes: T=" << std::setw(8) << std::fixed<< T + << " Size=" << std::setw(8) << usage.ru_maxrss/1024/1024 + << std::endl; +#endif + } + } + } + } + + return ret; + } + + void FieldGenerator::file_volumes(BasisClasses::BasisPtr basis, + CoefClasses::CoefsPtr coefs, + const std::string prefix, + const std::string outdir) + { + auto db = volumes(basis, coefs); + + if (myid==0) { + + int icnt = 0; + + for (auto & frame : db) { + + DataGrid datagrid(grid[0], grid[1], grid[2], + pmin[0], pmax[0], + pmin[1], pmax[1], + pmin[2], pmax[2]); + + std::vector tmp(grid[0]*grid[1]*grid[2]); + + for (auto & v : frame.second) { + + for (int i=0; i + FieldGenerator::histogram2d(PR::PRptr reader, std::vector ctr) + { + + std::map ret; + std::map fac; + + std::vector del(3, 0.0); + for (int k=0; k<3; k++) { + if (grid[k]>0) del[k] = (pmax[k] - pmin[k])/grid[k]; + } + + if (grid[0]>0 and grid[1]>0) { + ret["xy"] = Eigen::MatrixXf::Zero(grid[0], grid[1]); + fac["xy"] = 1.0/(grid[0]*grid[1]); + } + + if (grid[0]>0 and grid[2]>0) { + ret["xz"] = Eigen::MatrixXf::Zero(grid[0], grid[2]); + fac["xz"] = 1.0/(grid[0]*grid[2]); + } + + if (grid[1]>0 and grid[2]>0) { + ret["yz"] = Eigen::MatrixXf::Zero(grid[1], grid[2]); + fac["yz"] = 1.0/(grid[1]*grid[2]); + } + + std::vector pp(3); + std::vector bb(3); + + for (auto p=reader->firstParticle(); p!=0; p=reader->nextParticle()) { + + for (int k=0; k<3; k++) { + pp[k] = p->pos[k] - ctr[k]; + bb[k] = pp[k] >= pmin[k] and pp[k] < pmax[k] and del[k] > 0.0; + } + + // x-y + if (bb[0] and bb[1]) { + int indx1 = floor( (pp[0] - pmin[0])/del[0] ); + int indx2 = floor( (pp[1] - pmin[1])/del[1] ); + + if (indx1>=0 and indx1=0 and indx2mass * fac["xy"]; + } + + // x-z + if (bb[0] and bb[2]) { + int indx1 = floor( (pp[0] - pmin[0])/del[0] ); + int indx2 = floor( (pp[2] - pmin[2])/del[2] ); + + if (indx1>=0 and indx1=0 and indx2mass * fac["xz"]; + } + + // y-z + if (bb[1] and bb[2]) { + int indx1 = floor( (pp[1] - pmin[1])/del[1] ); + int indx2 = floor( (pp[2] - pmin[2])/del[2] ); + + if (indx1>=0 and indx1=0 and indx2mass * fac["yz"]; + } + + } + + if (use_mpi) { + for (auto & v : ret) { + if (myid==0) + MPI_Reduce(MPI_IN_PLACE, v.second.data(), v.second.size(), + MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD); + else + MPI_Reduce(v.second.data(), NULL, v.second.size(), + MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD); + } + } + + return ret; + } + // END histogram + + Eigen::VectorXf + FieldGenerator::histogram1d(PR::PRptr reader, double rmax, int nbins, + std::string proj, std::vector ctr) + { + const double pi = 3.14159265358979323846; + + Eigen::VectorXf ret = Eigen::VectorXf::Zero(nbins); + double del = rmax/nbins; + + enum class Projection {xy, xz, yz, r} type; + if (proj == "xy") type = Projection::xy; + else if (proj == "xz") type = Projection::xz; + else if (proj == "yz") type = Projection::yz; + else if (proj == "r" ) type = Projection::r; + else { + std::ostringstream sout; + sout << "FieldGenerator::histogram1d: error parsing projection <" << proj + << ">. Must be one of \"xy\", \"xz\", \"yz, \"r\"."; + std::runtime_error(sout.str()); + } + + // Make the histogram + // + for (auto p=reader->firstParticle(); p!=0; p=reader->nextParticle()) { + double rad = 0.0; + for (int k=0; k<3; k++) { + double pp = p->pos[k] - ctr[k]; + if (type == Projection::xy and (k==0 or k==1)) + rad += pp*pp; // Cylindrical radius x^2 + y^2 + else if (type == Projection::xz and (k==0 or k==2)) + rad += pp*pp; // Cylindrical radius x^2 + z^2 + else if (type == Projection::yz and (k==1 or k==2)) + rad += pp*pp; // Cylindrical radius y^2 + z^2 + else if (type == Projection::r) + rad += pp*pp; // Spherical radius + } + + int indx = floor(sqrt(rad)/del); + if (indx>=0 and indxmass; + } + + // Accumulate between MPI nodes; return value to root node + // + if (use_mpi) { + if (myid==0) + MPI_Reduce(MPI_IN_PLACE, ret.data(), ret.size(), + MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD); + else + MPI_Reduce(ret.data(), NULL, ret.size(), + MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD); + } + + // Inverse area or volume for density norm + // + for (int i=0; i +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace Utility +{ + + std::vector getDensityCenter(PR::PRptr reader, int nball, int Nsort) + { + std::vector ctr(3, 0.0); + + typedef point point3; + typedef kdtree tree3; + + std::vector points; + + double KDmass = 0.0; + for (auto p=reader->firstParticle(); p!=0; p=reader->nextParticle()) { + KDmass += p.mass; + points.push_back(point3({p.pos[0], p.pos[1], p.pos[2]}, p.mass)); + } + + std::vector dd; + int sz0 = points.size(); + + for (int n=0; n 0.0) { + + if (Nsort>0) { + + std::multimap stack, combo; + + for (int i=0; i(ret), 3.0); + if (volume>0.0) { + double density = std::get<1>(ret)/volume/KDmass; + stack[density] = i; + } + + if (stack.size()>Nsort) stack.erase(stack.begin()); + } + } + + combo = stack; + + for (int n=0; n dens; + std::vector indx; + for (auto v : stack) { + dens.push_back(v.first); + indx.push_back(v.second); + } + MPI_Bcast(&sz, 1, MPI_INT, n, MPI_COMM_WORLD); + MPI_Bcast(dens.data(), sz, MPI_DOUBLE, n, MPI_COMM_WORLD); + MPI_Bcast(indx.data(), sz, MPI_INT, n, MPI_COMM_WORLD); + } else { + int sz; + MPI_Bcast(&sz, 1, MPI_INT, n, MPI_COMM_WORLD); + std::vector dens(sz); + std::vector indx(sz); + MPI_Bcast(dens.data(), sz, MPI_DOUBLE, n, MPI_COMM_WORLD); + MPI_Bcast(indx.data(), sz, MPI_INT, n, MPI_COMM_WORLD); + for (int i=0; i Nsort) break; + for (int k=0; k<3; k++) + ctr[k] += it->first * points[it->second].get(k); + mastot += it->first; + } + + } else { + + for (int i=0; i(ret), 3.0); + double density = 0.0; + if (volume>0.0 and KDmass>0.0) + density = std::get<1>(ret)/volume/KDmass; + for (int k=0; k<3; k++) ctr[k] += density * points[i].get(k); + mastot += density; + } + } + + MPI_Allreduce(MPI_IN_PLACE, com.data(), 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &mastot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + } + + if (mastot>0.0) { + for (int k=0; k<3; k++) ctr[k] /= mastot; + } + } + + if (myid==0) { + std::cout << std::endl << "Computed density weighted center: ["; + for (int k=0; k<3; k++) std::cout << std::setw(16) << ctr[k]; + std::cout << "]" << std::endl; + } + } + + std::vector getCenterOfMass(PR::PRptr reader) + { + std::vector ctr(3, 0.0); + double mastot = 0.0; + + for (auto p=reader->firstParticle(); p!=0; p=reader->nextParticle()) { + for (int k=0; k<3; k++) ctr[k] += p.mass * p.pos[k]; + mastot += i.mass; + } + + MPI_Allreduce(MPI_IN_PLACE, com.data(), 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &mastot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + if (mastot>0.0) { + for (int k=0; k<3; k++) ctr[k] /= mastot; + } + + if (myid==0) { + std::cout << std::endl << "Computed COM: ["; + for (int k=0; k<3; k++) std::cout << std::setw(16) << ctr[k]; + std::cout << "]" << std::endl; + } + + return ctr; + } + +} +// END namespace Utility + diff --git a/expui/SvdSignChoice.o b/expui/SvdSignChoice.o new file mode 100644 index 000000000..44e447e46 Binary files /dev/null and b/expui/SvdSignChoice.o differ diff --git a/expui/matrix.tau4i.mat b/expui/matrix.tau4i.mat new file mode 100644 index 000000000..f95d347ed --- /dev/null +++ b/expui/matrix.tau4i.mat @@ -0,0 +1,417 @@ + 20 20 + 0 -0.5 1.04845 1.0065 -0.293603 + 0.210526 -0.5 0.989101 0.925431 -0.349137 + 0.421053 -0.5 0.868552 0.804073 -0.328403 + 0.631579 -0.5 0.772574 0.734059 -0.24089 + 0.842105 -0.5 0.740213 0.723049 -0.158479 + 1.05263 -0.5 0.741271 0.734274 -0.101611 + 1.26316 -0.5 0.750944 0.748414 -0.0615812 + 1.47368 -0.5 0.762483 0.761943 -0.0286814 + 1.68421 -0.5 0.777291 0.777291 0.00102937 + 1.89474 -0.5 0.796112 0.79565 0.0270991 + 2.10526 -0.5 0.818715 0.81719 0.0499565 + 2.31579 -0.5 0.847043 0.844327 0.0677793 + 2.52632 -0.5 0.877499 0.874334 0.0744592 + 2.73684 -0.5 0.901225 0.898431 0.0709001 + 2.94737 -0.5 0.916546 0.914267 0.0645911 + 3.15789 -0.5 0.926986 0.925087 0.0593053 + 3.36842 -0.5 0.934846 0.933196 0.055514 + 3.57895 -0.5 0.941635 0.940134 0.0531443 + 3.78947 -0.5 0.948231 0.946816 0.0517943 + 4 -0.5 0.955629 0.954297 0.0504464 + 0 -0.368421 1.10999 1.06398 -0.316262 + 0.210526 -0.368421 1.05573 0.957865 -0.443911 + 0.421053 -0.368421 0.850864 0.744224 -0.412433 + 0.631579 -0.368421 0.710032 0.661672 -0.257558 + 0.842105 -0.368421 0.692458 0.676358 -0.148453 + 1.05263 -0.368421 0.711123 0.705511 -0.0891668 + 1.26316 -0.368421 0.728769 0.726922 -0.0518516 + 1.47368 -0.368421 0.742309 0.742049 -0.0196344 + 1.68421 -0.368421 0.759221 0.759119 0.0124693 + 1.89474 -0.368421 0.781031 0.780011 0.0399093 + 2.10526 -0.368421 0.805168 0.802493 0.0655856 + 2.31579 -0.368421 0.839836 0.835078 0.0892788 + 2.52632 -0.368421 0.880882 0.875781 0.0946591 + 2.73684 -0.368421 0.908104 0.904266 0.0834017 + 2.94737 -0.368421 0.921775 0.918947 0.0721539 + 3.15789 -0.368421 0.930755 0.928501 0.0647383 + 3.36842 -0.368421 0.937102 0.935198 0.059706 + 3.57895 -0.368421 0.943142 0.941414 0.0570749 + 3.78947 -0.368421 0.948732 0.947073 0.0560718 + 4 -0.368421 0.957295 0.955668 0.0558002 + 0 -0.236842 1.20803 1.16821 -0.307573 + 0.210526 -0.236842 1.22942 1.06175 -0.619803 + 0.421053 -0.236842 0.8145 0.604634 -0.545735 + 0.631579 -0.236842 0.606705 0.554184 -0.246925 + 0.842105 -0.236842 0.638026 0.626885 -0.118708 + 1.05263 -0.236842 0.684152 0.680681 -0.0688304 + 1.26316 -0.236842 0.709318 0.708161 -0.0404998 + 1.47368 -0.236842 0.721166 0.721076 -0.011399 + 1.68421 -0.236842 0.739547 0.739067 0.0266434 + 1.89474 -0.236842 0.766278 0.764429 0.0532088 + 2.10526 -0.236842 0.787806 0.783687 0.080453 + 2.31579 -0.236842 0.829771 0.821121 0.119499 + 2.52632 -0.236842 0.892014 0.883611 0.122151 + 2.73684 -0.236842 0.920919 0.916011 0.0949418 + 2.94737 -0.236842 0.928529 0.925242 0.078053 + 3.15789 -0.236842 0.935845 0.933238 0.069798 + 3.36842 -0.236842 0.938795 0.936718 0.0624067 + 3.57895 -0.236842 0.94574 0.943801 0.0605289 + 3.78947 -0.236842 0.947414 0.945542 0.059527 + 4 -0.236842 0.95945 0.957336 0.0636465 + 0 -0.105263 1.32885 1.31411 -0.197386 + 0.210526 -0.105263 1.81865 1.53347 -0.977723 + 0.421053 -0.105263 0.704159 0.191183 -0.677709 + 0.631579 -0.105263 0.451974 0.422546 -0.160421 + 0.842105 -0.105263 0.591103 0.587857 -0.0618585 + 1.05263 -0.105263 0.66871 0.667377 -0.0422054 + 1.26316 -0.105263 0.695658 0.695057 -0.0289007 + 1.47368 -0.105263 0.697803 0.697723 -0.0105522 + 1.68421 -0.105263 0.716015 0.714166 0.0514266 + 1.89474 -0.105263 0.754653 0.751807 0.0654745 + 2.10526 -0.105263 0.768974 0.764524 0.08261 + 2.31579 -0.105263 0.811918 0.793546 0.171741 + 2.52632 -0.105263 0.926983 0.912082 0.16554 + 2.73684 -0.105263 0.948033 0.942768 0.099779 + 2.94737 -0.105263 0.938514 0.934956 0.0816436 + 3.15789 -0.105263 0.94471 0.941158 0.0818549 + 3.36842 -0.105263 0.934037 0.932384 0.0555299 + 3.57895 -0.105263 0.959145 0.957186 0.0612788 + 3.78947 -0.105263 0.941905 0.940115 0.0580474 + 4 -0.105263 0.96517 0.961898 0.0794108 + 0 0.0263158 1.37346 1.37235 0.0552535 + 0.210526 0.0263158 2.85474 2.74129 0.796812 + 0.421053 0.0263158 0.537829 -0.292095 0.451598 + 0.631579 0.0263158 0.367386 0.363986 0.0498686 + 0.842105 0.0263158 0.577456 0.577223 0.016406 + 1.05263 0.0263158 0.668534 0.66812 0.0235236 + 1.26316 0.0263158 0.697952 0.697721 0.0179675 + 1.47368 0.0263158 0.669303 0.669015 0.0196218 + 1.68421 0.0263158 0.68636 0.679147 -0.0992482 + 1.89474 0.0263158 0.764113 0.752927 -0.130267 + 2.10526 0.0263158 0.784747 0.783795 -0.0386266 + 2.31579 0.0263158 0.80073 0.754235 -0.268882 + 2.52632 0.0263158 1.05242 1.00775 -0.303363 + 2.73684 0.0263158 1.01105 1.00897 -0.0648253 + 2.94737 0.0263158 0.965166 0.96299 -0.0647776 + 3.15789 0.0263158 0.983963 0.97173 -0.154675 + 3.36842 0.0263158 0.925798 0.925548 -0.0215074 + 3.57895 0.0263158 0.989043 0.988571 -0.0305624 + 3.78947 0.0263158 0.944771 0.943799 -0.0428474 + 4 0.0263158 0.985661 0.982535 -0.078439 + 0 0.157895 1.28183 1.25522 0.259832 + 0.210526 0.157895 1.48407 1.24485 0.807971 + 0.421053 0.157895 0.765176 0.413586 0.643771 + 0.631579 0.157895 0.51848 0.474866 0.208144 + 0.842105 0.157895 0.607412 0.601011 0.0879545 + 1.05263 0.157895 0.672736 0.670628 0.0532155 + 1.26316 0.157895 0.700218 0.699414 0.0335462 + 1.47368 0.157895 0.707604 0.707546 0.00906002 + 1.68421 0.157895 0.726221 0.725184 -0.038791 + 1.89474 0.157895 0.758707 0.756275 -0.0607019 + 2.10526 0.157895 0.775337 0.770632 -0.0852883 + 2.31579 0.157895 0.820213 0.807073 -0.146226 + 2.52632 0.157895 0.907409 0.895809 -0.144626 + 2.73684 0.157895 0.933575 0.92827 -0.099384 + 2.94737 0.157895 0.933257 0.929788 -0.0803902 + 3.15789 0.157895 0.940145 0.937197 -0.0743994 + 3.36842 0.157895 0.9376 0.935601 -0.061187 + 3.57895 0.157895 0.950518 0.948482 -0.0621793 + 3.78947 0.157895 0.944475 0.942573 -0.059916 + 4 0.157895 0.961719 0.95904 -0.0717274 + 0 0.289474 1.16406 1.11976 0.318088 + 0.210526 0.289474 1.13776 1.00439 0.534503 + 0.421053 0.289474 0.83306 0.677046 0.485384 + 0.631579 0.289474 0.654267 0.601717 0.256909 + 0.842105 0.289474 0.660036 0.646379 0.133575 + 1.05263 0.289474 0.69413 0.689742 0.0779244 + 1.26316 0.289474 0.716582 0.715159 0.0451367 + 1.47368 0.289474 0.729787 0.729645 0.0143704 + 1.68421 0.289474 0.747681 0.747401 -0.0204514 + 1.89474 0.289474 0.771998 0.77051 -0.0479042 + 2.10526 0.289474 0.795348 0.79181 -0.0749442 + 2.31579 0.289474 0.834405 0.827656 -0.105911 + 2.52632 0.289474 0.886129 0.879275 -0.109995 + 2.73684 0.289474 0.914879 0.910375 -0.0906672 + 2.94737 0.289474 0.925702 0.922579 -0.0759802 + 3.15789 0.289474 0.933603 0.931144 -0.0677236 + 3.36842 0.289474 0.938343 0.936316 -0.0616489 + 3.57895 0.289474 0.944383 0.942526 -0.0591966 + 3.78947 0.289474 0.948324 0.946527 -0.0583586 + 4 0.289474 0.958479 0.956601 -0.0599802 + 0 0.421053 1.08157 1.03655 0.308794 + 0.210526 0.421053 1.02201 0.940578 0.399782 + 0.421053 0.421053 0.859098 0.773352 0.374134 + 0.631579 0.421053 0.738791 0.694317 0.252459 + 0.842105 0.421053 0.712689 0.695786 0.154297 + 1.05263 0.421053 0.723123 0.71685 0.0950437 + 1.26316 0.421053 0.737466 0.735336 0.0560114 + 1.47368 0.421053 0.750467 0.750105 0.0233009 + 1.68421 0.421053 0.766598 0.76656 -0.00767411 + 1.89474 0.421053 0.787119 0.786354 -0.0346767 + 2.10526 0.421053 0.810924 0.808759 -0.0592167 + 2.31579 0.421053 0.842887 0.839089 -0.0799325 + 2.52632 0.421053 0.87893 0.874721 -0.0859123 + 2.73684 0.421053 0.90479 0.901389 -0.0783701 + 2.94737 0.421053 0.919462 0.91685 -0.069255 + 3.15789 0.421053 0.929114 0.927 -0.0626463 + 3.36842 0.421053 0.936181 0.934375 -0.0581314 + 3.57895 0.421053 0.942502 0.940863 -0.0555612 + 3.78947 0.421053 0.948641 0.94708 -0.0543931 + 4 0.421053 0.956597 0.955099 -0.0535021 + 0 0.552632 1.03148 0.992042 0.282487 + 0.210526 0.552632 0.974647 0.9197 0.322627 + 0.421053 0.552632 0.873675 0.81913 0.303866 + 0.631579 0.552632 0.790395 0.755563 0.232052 + 0.842105 0.552632 0.75659 0.739695 0.158997 + 1.05263 0.552632 0.753122 0.745814 0.104662 + 1.26316 0.552632 0.760012 0.757246 0.0647866 + 1.47368 0.552632 0.770345 0.769677 0.0320738 + 1.68421 0.552632 0.784198 0.784192 0.00308834 + 1.89474 0.552632 0.801947 0.801638 -0.022277 + 2.10526 0.552632 0.823496 0.822314 -0.0441065 + 2.31579 0.552632 0.84967 0.847503 -0.0606521 + 2.52632 0.552632 0.877253 0.874636 -0.0677099 + 2.73684 0.552632 0.899603 0.897173 -0.0660791 + 2.94737 0.552632 0.91498 0.91292 -0.061366 + 3.15789 0.552632 0.925779 0.924025 -0.0569614 + 3.36842 0.552632 0.934027 0.932485 -0.0536451 + 3.57895 0.552632 0.941092 0.939684 -0.0514507 + 3.78947 0.552632 0.947863 0.946541 -0.0500447 + 4 0.552632 0.955025 0.953789 -0.0485686 + 0 0.684211 1.00177 0.968749 0.255102 + 0.210526 0.684211 0.95337 0.913568 0.272595 + 0.421053 0.684211 0.884213 0.846097 0.256813 + 0.631579 0.684211 0.823893 0.796749 0.209742 + 0.842105 0.684211 0.791166 0.775759 0.155374 + 1.05263 0.684211 0.780873 0.773284 0.108603 + 1.26316 0.684211 0.782228 0.779002 0.0709703 + 1.47368 0.684211 0.789375 0.788382 0.0395856 + 1.68421 0.684211 0.800721 0.800625 0.0123679 + 1.89474 0.684211 0.815857 0.815781 -0.0111737 + 2.10526 0.684211 0.834471 0.833901 -0.0308353 + 2.31579 0.684211 0.856068 0.854863 -0.0454032 + 2.52632 0.684211 0.878269 0.876658 -0.0531807 + 2.73684 0.684211 0.897506 0.895831 -0.0548099 + 2.94737 0.684211 0.912299 0.910745 -0.0532298 + 3.15789 0.684211 0.923473 0.922072 -0.0508459 + 3.36842 0.684211 0.932312 0.93104 -0.0486897 + 3.57895 0.684211 0.939855 0.938678 -0.0470205 + 3.78947 0.684211 0.946838 0.945735 -0.0456795 + 4 0.684211 0.95364 0.952615 -0.0441857 + 0 0.815789 0.984111 0.956718 0.230578 + 0.210526 0.815789 0.943718 0.913429 0.237174 + 0.421053 0.815789 0.892923 0.864573 0.223213 + 0.631579 0.815789 0.847336 0.825862 0.189549 + 0.842105 0.815789 0.818231 0.804711 0.148125 + 1.05263 0.815789 0.805246 0.79786 0.108816 + 1.26316 0.815789 0.802988 0.799499 0.0747793 + 1.47368 0.815789 0.807275 0.805989 0.0455565 + 1.68421 0.815789 0.816163 0.815912 0.0202247 + 1.89474 0.815789 0.828782 0.82878 -0.00153554 + 2.10526 0.815789 0.844501 0.844276 -0.0195004 + 2.31579 0.815789 0.862427 0.861798 -0.0329443 + 2.52632 0.815789 0.880782 0.879818 -0.0411835 + 2.73684 0.815789 0.897391 0.896275 -0.0447449 + 2.94737 0.815789 0.911111 0.909982 -0.0453364 + 3.15789 0.815789 0.922111 0.921032 -0.0446119 + 3.36842 0.815789 0.931115 0.930097 -0.0435257 + 3.57895 0.815789 0.93885 0.93789 -0.0424469 + 3.78947 0.815789 0.945864 0.944959 -0.0413703 + 4 0.815789 0.952438 0.951594 -0.0400698 + 0 0.947368 0.973608 0.950815 0.209437 + 0.210526 0.947368 0.939796 0.915919 0.210495 + 0.421053 0.947368 0.900576 0.878559 0.197919 + 0.631579 0.947368 0.864873 0.847559 0.172191 + 0.842105 0.947368 0.839709 0.828007 0.139694 + 1.05263 0.947368 0.826184 0.81925 0.106814 + 1.26316 0.947368 0.821795 0.818208 0.0767032 + 1.47368 0.947368 0.823812 0.82229 0.050043 + 1.68421 0.947368 0.830483 0.830054 0.0267077 + 1.89474 0.947368 0.840772 0.840746 0.00664627 + 2.10526 0.947368 0.853867 0.853809 -0.00992452 + 2.31579 0.947368 0.868796 0.868502 -0.022585 + 2.52632 0.947368 0.884188 0.883641 -0.0310782 + 2.73684 0.947368 0.898575 0.897861 -0.03583 + 2.94737 0.947368 0.911069 0.91028 -0.0379195 + 3.15789 0.947368 0.921569 0.920765 -0.0384931 + 3.36842 0.947368 0.930435 0.929645 -0.0383426 + 3.57895 0.947368 0.938137 0.937373 -0.0378616 + 3.78947 0.947368 0.945064 0.944334 -0.037149 + 4 0.947368 0.951436 0.95075 -0.0361347 + 0 1.07895 0.967465 0.948357 0.19133 + 0.210526 1.07895 0.938887 0.919567 0.189486 + 0.421053 1.07895 0.907481 0.889837 0.178079 + 0.631579 1.07895 0.878702 0.864482 0.157444 + 0.842105 1.07895 0.857095 0.846988 0.131235 + 1.05263 1.07895 0.84405 0.837666 0.103617 + 1.26316 1.07895 0.838539 0.834974 0.0772496 + 1.47368 1.07895 0.838876 0.837185 0.0532392 + 1.68421 1.07895 0.843663 0.843058 0.0319245 + 1.89474 1.07895 0.851884 0.851778 0.013476 + 2.10526 1.07895 0.862683 0.86268 -0.00188111 + 2.31579 1.07895 0.87513 0.87502 -0.0138895 + 2.52632 1.07895 0.888134 0.887849 -0.0224687 + 2.73684 1.07895 0.900621 0.900187 -0.027956 + 2.94737 1.07895 0.911883 0.911354 -0.0310758 + 3.15789 1.07895 0.921707 0.921129 -0.0326384 + 3.36842 1.07895 0.930234 0.929639 -0.0332751 + 3.57895 1.07895 0.937742 0.937148 -0.0333588 + 3.78947 1.07895 0.944491 0.943912 -0.033045 + 4 1.07895 0.950649 0.950099 -0.032356 + 0 1.21053 0.964039 0.947882 0.175762 + 0.210526 1.21053 0.939624 0.923673 0.172397 + 0.421053 1.21053 0.913782 0.899304 0.162016 + 0.631579 1.21053 0.89004 0.878167 0.144888 + 0.842105 1.21053 0.871453 0.862694 0.123241 + 1.05263 1.21053 0.859313 0.853491 0.0998562 + 1.26316 1.21053 0.853319 0.849852 0.0768428 + 1.47368 1.21053 0.852468 0.850668 0.0553775 + 1.68421 1.21053 0.855717 0.854958 0.0360256 + 1.89474 1.21053 0.862164 0.861952 0.0191009 + 2.10526 1.21053 0.870993 0.87098 0.00484675 + 2.31579 1.21053 0.881367 0.881342 -0.00655454 + 2.52632 1.21053 0.892392 0.892265 -0.015087 + 2.73684 1.21053 0.90324 0.902995 -0.0210102 + 2.94737 1.21053 0.913324 0.912987 -0.0248306 + 3.15789 1.21053 0.922393 0.921994 -0.0271331 + 3.36842 1.21053 0.930457 0.930023 -0.0284136 + 3.57895 1.21053 0.937658 0.937209 -0.0290052 + 3.78947 1.21053 0.94416 0.943712 -0.0290861 + 4 1.21053 0.950082 0.949647 -0.0287308 + 0 1.34211 0.962344 0.948562 0.16228 + 0.210526 1.34211 0.941268 0.927886 0.158151 + 0.421053 1.34211 0.919558 0.907457 0.148687 + 0.631579 1.34211 0.899601 0.889547 0.13412 + 0.842105 1.34211 0.883523 0.875889 0.115889 + 1.05263 1.34211 0.872417 0.867129 0.095904 + 1.26316 1.34211 0.866321 0.862998 0.0758039 + 1.47368 1.34211 0.864665 0.862805 0.0566804 + 1.68421 1.34211 0.866692 0.865806 0.0391764 + 1.89474 1.34211 0.871652 0.87133 0.0236808 + 2.10526 1.34211 0.878815 0.878753 0.0104528 + 2.31579 1.34211 0.887447 0.887447 -0.000353893 + 2.52632 1.34211 0.896808 0.896765 -0.00873672 + 2.73684 1.34211 0.906228 0.906106 -0.0148892 + 2.94737 1.34211 0.915217 0.915016 -0.0191735 + 3.15789 1.34211 0.92351 0.923247 -0.022019 + 3.36842 1.34211 0.931042 0.930737 -0.0238141 + 3.57895 1.34211 0.937865 0.937536 -0.0248462 + 3.78947 1.34211 0.944068 0.943729 -0.0252959 + 4 1.34211 0.94973 0.949394 -0.0252639 + 0 1.47368 0.961762 0.949912 0.150514 + 0.210526 1.47368 0.943401 0.932027 0.146051 + 0.421053 1.47368 0.924863 0.914598 0.137413 + 0.631579 1.47368 0.907829 0.89921 0.124804 + 0.842105 1.47368 0.893824 0.887128 0.109206 + 1.05263 1.47368 0.88374 0.878941 0.0919693 + 1.26316 1.47368 0.87776 0.874605 0.0743628 + 1.47368 1.47368 0.875578 0.873698 0.0573397 + 1.68421 1.47368 0.876657 0.875672 0.0415376 + 1.89474 1.47368 0.880391 0.879966 0.0273706 + 2.10526 1.47368 0.886157 0.886028 0.0151062 + 2.31579 1.47368 0.893323 0.89331 0.00489047 + 2.52632 1.47368 0.90127 0.901264 -0.00326591 + 2.73684 1.47368 0.909445 0.909396 -0.00950104 + 2.94737 1.47368 0.917427 0.917319 -0.0140762 + 3.15789 1.47368 0.924957 0.924795 -0.0173103 + 3.36842 1.47368 0.931928 0.931724 -0.0195074 + 3.57895 1.47368 0.93833 0.938097 -0.0209111 + 3.78947 1.47368 0.944201 0.943952 -0.0216921 + 4 1.47368 0.949583 0.949329 -0.0219616 + 0 1.60526 0.961904 0.951637 0.140167 + 0.210526 1.60526 0.945778 0.936004 0.135621 + 0.421053 1.60526 0.92974 0.920925 0.127728 + 0.631579 1.60526 0.915016 0.907547 0.116673 + 0.842105 1.60526 0.902728 0.896814 0.103157 + 1.05263 1.60526 0.893593 0.889234 0.0881633 + 1.26316 1.60526 0.887846 0.884867 0.072679 + 1.47368 1.60526 0.885337 0.883467 0.0575121 + 1.68421 1.60526 0.885689 0.884633 0.0432539 + 1.89474 1.60526 0.888427 0.88791 0.030312 + 2.10526 1.60526 0.893028 0.892827 0.0189532 + 2.31579 1.60526 0.898959 0.898911 0.00932338 + 2.52632 1.60526 0.9057 0.905699 0.00144852 + 2.73684 1.60526 0.912786 0.912774 -0.00476394 + 2.94737 1.60526 0.919851 0.919802 -0.00950196 + 3.15789 1.60526 0.926652 0.926561 -0.0130036 + 3.36842 1.60526 0.933057 0.932928 -0.0155072 + 3.57895 1.60526 0.939019 0.938862 -0.0172164 + 3.78947 1.60526 0.944538 0.944361 -0.0182865 + 4 1.60526 0.949626 0.94944 -0.0188292 + 0 1.73684 0.962514 0.953556 0.131009 + 0.210526 1.73684 0.948249 0.93977 0.126523 + 0.421053 1.73684 0.934227 0.926577 0.119304 + 0.631579 1.73684 0.921363 0.914831 0.109517 + 0.842105 1.73684 0.910503 0.905248 0.0976831 + 1.05263 1.73684 0.902226 0.898256 0.0845401 + 1.26316 1.73684 0.896768 0.893964 0.0708608 + 1.47368 1.73684 0.89407 0.89223 0.0573213 + 1.68421 1.73684 0.893874 0.892768 0.0444498 + 1.89474 1.73684 0.895807 0.895213 0.0326298 + 2.10526 1.73684 0.899443 0.899171 0.0221193 + 2.31579 1.73684 0.904334 0.90424 0.0130652 + 2.52632 1.73684 0.910042 0.910026 0.00550909 + 2.73684 1.73684 0.916175 0.916175 -0.000605093 + 2.94737 1.73684 0.922412 0.922396 -0.00541056 + 3.15789 1.73684 0.928527 0.928483 -0.00908502 + 3.36842 1.73684 0.934379 0.934304 -0.0118151 + 3.57895 1.73684 0.939897 0.939796 -0.0137693 + 3.78947 1.73684 0.945053 0.944933 -0.0150855 + 4 1.73684 0.949841 0.949709 -0.0158702 + 0 1.86842 0.963422 0.955557 0.122851 + 0.210526 1.86842 0.950725 0.94331 0.11851 + 0.421053 1.86842 0.938354 0.931658 0.1119 + 0.631579 1.86842 0.927014 0.921255 0.10317 + 0.842105 1.86842 0.917353 0.912654 0.0927222 + 1.05263 1.86842 0.909837 0.906213 0.0811215 + 1.26316 1.86842 0.904688 0.902054 0.068981 + 1.47368 1.86842 0.901899 0.900104 0.0568628 + 1.68421 1.86842 0.901292 0.900156 0.0452294 + 1.89474 1.86842 0.902581 0.901924 0.034431 + 2.10526 1.86842 0.905419 0.905082 0.0247117 + 2.31579 1.86842 0.909434 0.90929 0.0162175 + 2.52632 1.86842 0.914258 0.914214 0.00900281 + 2.73684 1.86842 0.919556 0.919551 0.0030405 + 2.94737 1.86842 0.925049 0.925047 -0.0017613 + 3.15789 1.86842 0.930527 0.930511 -0.00553443 + 3.36842 1.86842 0.935849 0.935811 -0.00842453 + 3.57895 1.86842 0.94093 0.940871 -0.0105701 + 3.78947 1.86842 0.945724 0.945647 -0.012091 + 4 1.86842 0.950209 0.950119 -0.0130856 + 0 2 0.964515 0.957569 0.115546 + 0.210526 2 0.95315 0.946618 0.111395 + 0.421053 2 0.942154 0.936247 0.105336 + 0.631579 2 0.932078 0.926964 0.097503 + 0.842105 2 0.92343 0.919207 0.0882133 + 1.05263 2 0.916587 0.913269 0.0779111 + 1.26316 2 0.911747 0.909275 0.0670875 + 1.47368 2 0.908932 0.907193 0.0562099 + 1.68421 2 0.908023 0.906873 0.045678 + 1.89474 2 0.908798 0.908092 0.0358062 + 2.10526 2 0.910978 0.910583 0.0268214 + 2.31579 2 0.914256 0.914061 0.0188661 + 2.52632 2 0.918321 0.918243 0.0120043 + 2.73684 2 0.922889 0.922868 0.00623095 + 2.94737 2 0.927715 0.927714 0.00148556 + 3.15789 2 0.932608 0.932605 -0.00232839 + 3.36842 2 0.937429 0.937414 -0.00532388 + 3.57895 2 0.942087 0.942056 -0.00761388 + 3.78947 2 0.946525 0.94648 -0.0093012 + 4 2 0.95071 0.950652 -0.0104747 +--- +nptsE: 6 +nptsK: 5 +M: 2 +mmax: 4 +l1max: 10 +nmax: 40 +model: Zang +Itype: 1 +A: 1 +nu: 4 +mu: 5 +sigma: 0.37796447301000002 +scl: 1 +... + diff --git a/expui/patch.txt b/expui/patch.txt new file mode 100644 index 000000000..e397e4987 --- /dev/null +++ b/expui/patch.txt @@ -0,0 +1,120 @@ +diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc +index 25f20ab1..ab9945ca 100644 +--- a/expui/BiorthBasis.cc ++++ b/expui/BiorthBasis.cc +@@ -3885,11 +3885,29 @@ namespace BasisClasses + + // Sanity check + // ++ if (tfinal == tinit) { ++ throw std::runtime_error ++ ("BasisClasses::IntegrateOrbits: tinit cannot be equal to tfinal"); ++ } ++ ++ if (h < 0.0 and tfinal > tinit) { ++ std::ostringstream sout; ++ throw std::runtime_error ++ ("BasisClasses::IntegrateOrbits: tfinal must be smaller than tinit " ++ "when step size is negative"); ++ } ++ ++ if (h > 0.0 and tfinal < tinit) { ++ throw std::runtime_error ++ ("BasisClasses::IntegrateOrbits: tfinal must be larger than " ++ "tinit when step size is positive"); ++ } ++ + if ( (tfinal - tinit)/h > + static_cast(std::numeric_limits::max()) ) + { +- std::cout << "BasicFactor::IntegrateOrbits: step size is too small or " +- << "time interval is too large.\n"; ++ std::cout << "BasisClasses::IntegrateOrbits: step size is too small or " ++ << "time interval is too large." << std::endl; + // Return empty data + // + return {Eigen::VectorXd(), Eigen::Tensor()}; +@@ -3897,16 +3915,22 @@ namespace BasisClasses + + // Number of steps + // +- int numT = floor( (tfinal - tinit)/h ); ++ int numT = std::ceil( (tfinal - tinit)/h ); ++ ++ // Want both end points in the output at minimum ++ // ++ numT = std::max(2, numT); ++ nout = std::max(2, nout); + +- // Number of output steps. Will attempt to find the best stride... ++ // Number of output steps. Compute a stride>1 if nout(numT)/static_cast(nout)); +- if (stride>1) numT = nout * stride; ++ if (stride>1) numT = (nout-1) * stride + 1; ++ else nout = numT; + + // Compute the interval-matching step + // +- h = (tfinal - tinit)/(numT - 1); ++ h = (tfinal - tinit)/(numT-1); + + // Return data + // +@@ -3916,10 +3940,10 @@ namespace BasisClasses + ret.resize(rows, 6, nout); + } + catch (const std::bad_alloc& e) { +- std::cout << "BasicFactor::IntegrateOrbits: memory allocation failed: " ++ std::cout << "BasisClasses::IntegrateOrbits: memory allocation failed: " + << e.what() << std::endl + << "Your requested number of orbits and time steps requires " +- << floor(8.0*rows*6*nout/1e9)+1 << " GB free memory" ++ << floor(4.0*rows*6*nout/stride/1e9)+1 << " GB free memory" + << std::endl; + + // Return empty data +@@ -3931,18 +3955,24 @@ namespace BasisClasses + // + Eigen::VectorXd times(nout); + +- // Do the work ++ // Assign the initial point + // + times(0) = tinit; + for (int n=0; n +//
  • Savitzky-Golay smoothing. +//
  • computing a numerical derivative based of Savitzky-Golay smoothing. +//
  • required linear algebra support for SG smoothing using STL based +// vector/matrix classes +// +// +// \brief Linear Algebra "Toolkit". +// + +// system headers +#include +#include // for size_t +#include // for fabs +#include + +#include "sgsmooth.h" + +//! default convergence +static const double TINY_FLOAT = 1.0e-300; + +//! comfortable array of doubles +typedef std::vector float_vect; +//! comfortable array of ints; +typedef std::vector int_vect; + +// savitzky golay smoothing. +static float_vect sg_smooth(const float_vect &v, const int w, const int deg); +//! numerical derivative based on savitzky golay smoothing. +static float_vect sg_derivative(const float_vect &v, const int w, + const int deg, const double h=1.0); + +// c-callable interface to for tcl plugin code +extern "C" double *calc_sgsmooth(const int ndat, double *input, + const int window, const int order) +{ + float_vect in(ndat); + int i; + +#if defined(_OPENMP) +#pragma omp parallel for private(i) schedule(static) +#endif + for (i=0; i < ndat; ++i) { + in[i] = input[i]; + } + + float_vect out = sg_smooth(in, window, order); + +#if defined(_OPENMP) +#pragma omp parallel for private(i) schedule(static) +#endif + for (i=0; i < ndat; ++i) { + input[i] = out[i]; + } + + return input; +} + + +extern "C" double *calc_sgsderiv(const int ndat, double *input, + const int window, const int order, + const double delta) +{ + float_vect in(ndat); + int i; + +#if defined(_OPENMP) +#pragma omp parallel for private(i) schedule(static) +#endif + for (i=0; i < ndat; ++i) { + in[i] = input[i]; + } + + float_vect out = sg_derivative(in, window, order, delta); + +#if defined(_OPENMP) +#pragma omp parallel for private(i) schedule(static) +#endif + for (i=0; i < ndat; ++i) { + input[i] = out[i]; + } + + return input; +} + +/*! matrix class. + * + * This is a matrix class derived from a vector of float_vects. Note that + * the matrix elements indexed [row][column] with indices starting at 0 (c + * style). Also note that because of its design looping through rows should + * be faster than looping through columns. + * + * \brief two dimensional floating point array + */ +class float_mat : public std::vector { +private: + //! disable the default constructor + explicit float_mat() {}; + //! disable assignment operator until it is implemented. + float_mat &operator =(const float_mat &) { return *this; }; +public: + //! constructor with sizes + float_mat(const size_t rows, const size_t cols, const double def=0.0); + //! copy constructor for matrix + float_mat(const float_mat &m); + //! copy constructor for vector + float_mat(const float_vect &v); + + //! use default destructor + // ~float_mat() {}; + + //! get size + size_t nr_rows(void) const { return size(); }; + //! get size + size_t nr_cols(void) const { return front().size(); }; +}; + + + +// constructor with sizes +float_mat::float_mat(const size_t rows,const size_t cols,const double defval) + : std::vector(rows) { + int i; + for (i = 0; i < rows; ++i) { + (*this)[i].resize(cols, defval); + } + if ((rows < 1) || (cols < 1)) { + char buffer[1024]; + + sprintf(buffer, "cannot build matrix with %d rows and %d columns\n", + rows, cols); + sgs_error(buffer); + } +} + +// copy constructor for matrix +float_mat::float_mat(const float_mat &m) : std::vector(m.size()) { + + float_mat::iterator inew = begin(); + float_mat::const_iterator iold = m.begin(); + for (/* empty */; iold < m.end(); ++inew, ++iold) { + const size_t oldsz = iold->size(); + inew->resize(oldsz); + const float_vect oldvec(*iold); + *inew = oldvec; + } +} + +// copy constructor for vector +float_mat::float_mat(const float_vect &v) + : std::vector(1) { + + const size_t oldsz = v.size(); + front().resize(oldsz); + front() = v; +} + +////////////////////// +// Helper functions // +////////////////////// + +//! permute() orders the rows of A to match the integers in the index array. +void permute(float_mat &A, int_vect &idx) +{ + int_vect i(idx.size()); + int j,k; + + for (j = 0; j < A.nr_rows(); ++j) { + i[j] = j; + } + + // loop over permuted indices + for (j = 0; j < A.nr_rows(); ++j) { + if (i[j] != idx[j]) { + + // search only the remaining indices + for (k = j+1; k < A.nr_rows(); ++k) { + if (i[k] ==idx[j]) { + std::swap(A[j],A[k]); // swap the rows and + i[k] = i[j]; // the elements of + i[j] = idx[j]; // the ordered index. + break; // next j + } + } + } + } +} + +/*! \brief Implicit partial pivoting. + * + * The function looks for pivot element only in rows below the current + * element, A[idx[row]][column], then swaps that row with the current one in + * the index map. The algorithm is for implicit pivoting (i.e., the pivot is + * chosen as if the max coefficient in each row is set to 1) based on the + * scaling information in the vector scale. The map of swapped indices is + * recorded in swp. The return value is +1 or -1 depending on whether the + * number of row swaps was even or odd respectively. */ +static int partial_pivot(float_mat &A, const size_t row, const size_t col, + float_vect &scale, int_vect &idx, double tol) +{ + if (tol <= 0.0) + tol = TINY_FLOAT; + + int swapNum = 1; + + // default pivot is the current position, [row,col] + int pivot = row; + double piv_elem = fabs(A[idx[row]][col]) * scale[idx[row]]; + + // loop over possible pivots below current + int j; + for (j = row + 1; j < A.nr_rows(); ++j) { + + const double tmp = fabs(A[idx[j]][col]) * scale[idx[j]]; + + // if this elem is larger, then it becomes the pivot + if (tmp > piv_elem) { + pivot = j; + piv_elem = tmp; + } + } + +#if 0 + if(piv_elem < tol) { + sgs_error("partial_pivot(): Zero pivot encountered.\n") +#endif + + if(pivot > row) { // bring the pivot to the diagonal + j = idx[row]; // reorder swap array + idx[row] = idx[pivot]; + idx[pivot] = j; + swapNum = -swapNum; // keeping track of odd or even swap + } + return swapNum; +} + +/*! \brief Perform backward substitution. + * + * Solves the system of equations A*b=a, ASSUMING that A is upper + * triangular. If diag==1, then the diagonal elements are additionally + * assumed to be 1. Note that the lower triangular elements are never + * checked, so this function is valid to use after a LU-decomposition in + * place. A is not modified, and the solution, b, is returned in a. */ +static void lu_backsubst(float_mat &A, float_mat &a, bool diag=false) +{ + int r,c,k; + + for (r = (A.nr_rows() - 1); r >= 0; --r) { + for (c = (A.nr_cols() - 1); c > r; --c) { + for (k = 0; k < A.nr_cols(); ++k) { + a[r][k] -= A[r][c] * a[c][k]; + } + } + if(!diag) { + for (k = 0; k < A.nr_cols(); ++k) { + a[r][k] /= A[r][r]; + } + } + } +} + +/*! \brief Perform forward substitution. + * + * Solves the system of equations A*b=a, ASSUMING that A is lower + * triangular. If diag==1, then the diagonal elements are additionally + * assumed to be 1. Note that the upper triangular elements are never + * checked, so this function is valid to use after a LU-decomposition in + * place. A is not modified, and the solution, b, is returned in a. */ +static void lu_forwsubst(float_mat &A, float_mat &a, bool diag=true) +{ + int r,k,c; + for (r = 0;r < A.nr_rows(); ++r) { + for(c = 0; c < r; ++c) { + for (k = 0; k < A.nr_cols(); ++k) { + a[r][k] -= A[r][c] * a[c][k]; + } + } + if(!diag) { + for (k = 0; k < A.nr_cols(); ++k) { + a[r][k] /= A[r][r]; + } + } + } +} + +/*! \brief Performs LU factorization in place. + * + * This is Crout's algorithm (cf., Num. Rec. in C, Section 2.3). The map of + * swapped indeces is recorded in idx. The return value is +1 or -1 + * depending on whether the number of row swaps was even or odd + * respectively. idx must be preinitialized to a valid set of indices + * (e.g., {1,2, ... ,A.nr_rows()}). */ +static int lu_factorize(float_mat &A, int_vect &idx, double tol=TINY_FLOAT) +{ + if ( tol <= 0.0) + tol = TINY_FLOAT; + + if ((A.nr_rows() == 0) || (A.nr_rows() != A.nr_cols())) { + sgs_error("lu_factorize(): cannot handle empty " + "or nonsquare matrices.\n"); + + return 0; + } + + float_vect scale(A.nr_rows()); // implicit pivot scaling + int i,j; + for (i = 0; i < A.nr_rows(); ++i) { + double maxval = 0.0; + for (j = 0; j < A.nr_cols(); ++j) { + if (fabs(A[i][j]) > maxval) + maxval = fabs(A[i][j]); + } + if (maxval == 0.0) { + sgs_error("lu_factorize(): zero pivot found.\n"); + return 0; + } + scale[i] = 1.0 / maxval; + } + + int swapNum = 1; + int c,r; + for (c = 0; c < A.nr_cols() ; ++c) { // loop over columns + swapNum *= partial_pivot(A, c, c, scale, idx, tol); // bring pivot to diagonal + for(r = 0; r < A.nr_rows(); ++r) { // loop over rows + int lim = (r < c) ? r : c; + for (j = 0; j < lim; ++j) { + A[idx[r]][c] -= A[idx[r]][j] * A[idx[j]][c]; + } + if (r > c) + A[idx[r]][c] /= A[idx[c]][c]; + } + } + permute(A,idx); + return swapNum; +} + +/*! \brief Solve a system of linear equations. + * Solves the inhomogeneous matrix problem with lu-decomposition. Note that + * inversion may be accomplished by setting a to the identity_matrix. */ +static float_mat lin_solve(const float_mat &A, const float_mat &a, + double tol=TINY_FLOAT) +{ + float_mat B(A); + float_mat b(a); + int_vect idx(B.nr_rows()); + int j; + + for (j = 0; j < B.nr_rows(); ++j) { + idx[j] = j; // init row swap label array + } + lu_factorize(B,idx,tol); // get the lu-decomp. + permute(b,idx); // sort the inhomogeneity to match the lu-decomp + lu_forwsubst(B,b); // solve the forward problem + lu_backsubst(B,b); // solve the backward problem + return b; +} + +/////////////////////// +// related functions // +/////////////////////// + +//! Returns the inverse of a matrix using LU-decomposition. +static float_mat invert(const float_mat &A) +{ + const int n = A.size(); + float_mat E(n, n, 0.0); + float_mat B(A); + int i; + + for (i = 0; i < n; ++i) { + E[i][i] = 1.0; + } + + return lin_solve(B, E); +} + +//! returns the transposed matrix. +static float_mat transpose(const float_mat &a) +{ + float_mat res(a.nr_cols(), a.nr_rows()); + int i,j; + + for (i = 0; i < a.nr_rows(); ++i) { + for (j = 0; j < a.nr_cols(); ++j) { + res[j][i] = a[i][j]; + } + } + return res; +} + +//! matrix multiplication. +float_mat operator *(const float_mat &a, const float_mat &b) +{ + float_mat res(a.nr_rows(), b.nr_cols()); + if (a.nr_cols() != b.nr_rows()) { + sgs_error("incompatible matrices in multiplication\n"); + return res; + } + + int i,j,k; + + for (i = 0; i < a.nr_rows(); ++i) { + for (j = 0; j < b.nr_cols(); ++j) { + double sum(0.0); + for (k = 0; k < a.nr_cols(); ++k) { + sum += a[i][k] * b[k][j]; + } + res[i][j] = sum; + } + } + return res; +} + + +//! calculate savitzky golay coefficients. +static float_vect sg_coeff(const float_vect &b, const size_t deg) +{ + const size_t rows(b.size()); + const size_t cols(deg + 1); + float_mat A(rows, cols); + float_vect res(rows); + + // generate input matrix for least squares fit + int i,j; + for (i = 0; i < rows; ++i) { + for (j = 0; j < cols; ++j) { + A[i][j] = pow(double(i), double(j)); + } + } + + float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b))); + + for (i = 0; i < b.size(); ++i) { + res[i] = c[0][0]; + for (j = 1; j <= deg; ++j) { + res[i] += c[j][0] * pow(double(i), double(j)); + } + } + return res; +} + +/*! \brief savitzky golay smoothing. + * + * This method means fitting a polynome of degree 'deg' to a sliding window + * of width 2w+1 throughout the data. The needed coefficients are + * generated dynamically by doing a least squares fit on a "symmetric" unit + * vector of size 2w+1, e.g. for w=2 b=(0,0,1,0,0). evaluating the polynome + * yields the sg-coefficients. at the border non symmectric vectors b are + * used. */ +float_vect sg_smooth(const float_vect &v, const int width, const int deg) +{ + float_vect res(v.size(), 0.0); + if ((width < 1) || (deg < 1) || (v.size() < (2 * width + 2))) { + sgs_error("sgsmooth: parameter error.\n"); + return res; + } + + const int window = 2 * width + 1; + const int endidx = v.size() - 1; + + // handle border cases first because we need different coefficients + int i,j; +#if defined(_OPENMP) +#pragma omp parallel for private(i,j) schedule(static) +#endif + for (i = 0; i < width; ++i) { + float_vect b1(window, 0.0); + b1[i] = 1.0; + + const float_vect c1(sg_coeff(b1, deg)); + for (j = 0; j < window; ++j) { + res[i] += c1[j] * v[j]; + res[endidx - i] += c1[j] * v[endidx - j]; + } + } + + // now loop over rest of data. reusing the "symmetric" coefficients. + float_vect b2(window, 0.0); + b2[width] = 1.0; + const float_vect c2(sg_coeff(b2, deg)); + +#if defined(_OPENMP) +#pragma omp parallel for private(i,j) schedule(static) +#endif + for (i = 0; i <= (v.size() - window); ++i) { + for (j = 0; j < window; ++j) { + res[i + width] += c2[j] * v[i + j]; + } + } + return res; +} + +/*! least squares fit a polynome of degree 'deg' to data in 'b'. + * then calculate the first derivative and return it. */ +static float_vect lsqr_fprime(const float_vect &b, const int deg) +{ + const int rows(b.size()); + const int cols(deg + 1); + float_mat A(rows, cols); + float_vect res(rows); + + // generate input matrix for least squares fit + int i,j; + for (i = 0; i < rows; ++i) { + for (j = 0; j < cols; ++j) { + A[i][j] = pow(double(i), double(j)); + } + } + + float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b))); + + for (i = 0; i < b.size(); ++i) { + res[i] = c[1][0]; + for (j = 1; j < deg; ++j) { + res[i] += c[j + 1][0] * double(j+1) + * pow(double(i), double(j)); + } + } + return res; +} + +/*! \brief savitzky golay smoothed numerical derivative. + * + * This method means fitting a polynome of degree 'deg' to a sliding window + * of width 2w+1 throughout the data. + * + * In contrast to the sg_smooth function we do a brute force attempt by + * always fitting the data to a polynome of degree 'deg' and using the + * result. */ +float_vect sg_derivative(const float_vect &v, const int width, + const int deg, const double h) +{ + float_vect res(v.size(), 0.0); + if ((width < 1) || (deg < 1) || (v.size() < (2 * width + 2))) { + sgs_error("sgsderiv: parameter error.\n"); + return res; + } + + const int window = 2 * width + 1; + + // handle border cases first because we do not repeat the fit + // lower part + float_vect b(window, 0.0); + int i,j; + + for (i = 0; i < window; ++i) { + b[i] = v[i] / h; + } + const float_vect c(lsqr_fprime(b, deg)); + for (j = 0; j <= width; ++j) { + res[j] = c[j]; + } + // upper part. direction of fit is reversed + for (i = 0; i < window; ++i) { + b[i] = v[v.size() - 1 - i] / h; + } + const float_vect d(lsqr_fprime(b, deg)); + for (i = 0; i <= width; ++i) { + res[v.size() - 1 - i] = -d[i]; + } + + // now loop over rest of data. wasting a lot of least squares calcs + // since we only use the middle value. +#if defined(_OPENMP) +#pragma omp parallel for private(i,j) schedule(static) +#endif + for (i = 1; i < (v.size() - window); ++i) { + for (j = 0; j < window; ++j) { + b[j] = v[i + j] / h; + } + res[i + width] = lsqr_fprime(b, deg)[width]; + } + return res; +} + +// Local Variables: +// mode: c++ +// c-basic-offset: 4 +// fill-column: 76 +// indent-tabs-mode: nil +// End: + + diff --git a/expui/store.resize b/expui/store.resize new file mode 100644 index 000000000..fa48d37da --- /dev/null +++ b/expui/store.resize @@ -0,0 +1,45 @@ + cf->lmax = lmax; + cf->nmax = nmax; + cf->scale = scale; + cf->time = time; + cf->normed = true; + cf->store.resize((lmax+1)*(lmax+2)/2, nmax); + cf->coefs = std::make_shared + (cf->store.data(), ldim, nmax); + (*cf->coefs)(L0, n) = {expcoef(L1, n), 0.0}; + (*cf->coefs)(L0, n) = {expcoef(L1, n), expcoef(L1+1, n)}; + expcoef(L1, n) = (*cf->coefs)(L0, n).real(); + expcoef(L1, n) = (*cf->coefs)(L0, n).real(); + expcoef(L1+1, n) = (*cf->coefs)(L0, n).imag(); + if (cf->ctr.size()) + coefctr = cf->ctr; + cf->mmax = mmax; + cf->nmax = nmax; + cf->time = time; + cf->store.resize((mmax+1)*nmax); + cf->coefs = std::make_shared + (cf->store.data(), mmax+1, nmax); + (*cf->coefs)(m, n) = {cos1(n), sin1(n)}; + sl->set_coefs(m, (*cf->coefs).row(m).real(), (*cf->coefs).row(m).imag(), m==0); + if (cf->ctr.size()) + coefctr = cf->ctr; + cf->mmax = mmax; + cf->nmax = nmax; + cf->time = time; + cf->store((2*mmax+1)*nmax); + cf->coefs = std::make_shared + (cf->store.data(), 2*mmax+1, nmax); + (*cf->coefs)(m, n) = {expcoef(m0, n), 0.0}; + (*cf->coefs)(m, n) = {expcoef(m0, n), expcoef(m0+1, n)}; + int rows = cf->rows(); + int cols = cf->cols(); + auto & cc = *cf->coefs; + if (cf->ctr.size()) + coefctr = cf->ctr; + cf->nmaxx = nmaxx; + cf->nmaxy = nmaxy; + cf->nmaxz = nmaxz; + cf->time = time; + cf->allocate(); + *cf->coefs = expcoef; + expcoef = *cf->coefs; diff --git a/expui/svd_test b/expui/svd_test new file mode 100755 index 000000000..56bc392dd Binary files /dev/null and b/expui/svd_test differ diff --git a/expui/svd_test.cc b/expui/svd_test.cc new file mode 100644 index 000000000..63c22c301 --- /dev/null +++ b/expui/svd_test.cc @@ -0,0 +1,112 @@ +#include +#include + +namespace MSSA { + void SvdSignChoice + (const Eigen::MatrixXd& X, + Eigen::MatrixXd& U, const Eigen::VectorXd& S, Eigen::MatrixXd& V, + bool sample); +} + +int main() +{ + srand((unsigned int) time(0)); + + // Instantiate a matrix filled with random numbers + const int NROWS = 32; // N + const int MCOLS = 32; // M + + Eigen::MatrixXd A = Eigen::MatrixXd::Random(NROWS, MCOLS); + std::cout << "Here is the data matrix A:" << std::endl << A + << std::endl << std::endl; + + // Carry out the SVD decomposition of the matrix, A such that, + // A = U S V^T where A is NxM, and + // U is NxM, S=diag(w_{0}, w_{1} ...) is MxM, and V is MxM. + // U is column-orthogonal + // D is the diagonal matrix with the singular values + // V is orthogonal + + Eigen::JacobiSVD + svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV); + + std::cout << "Its singular values are:" << std::endl + << svd.singularValues() << std::endl << std::endl; + + std::cout << "Its left singular vectors are the columns of the thin U matrix:" + << std::endl << svd.matrixU() << std::endl << std::endl; + + std::cout << "Its right singular vectors are the columns of the thin V matrix:" + << std::endl << svd.matrixV() << std::endl << std::endl; + + // Test that this SVD is doing what we expect. + // + Eigen::MatrixXd S = svd.singularValues().asDiagonal(); + + std::cout << std::endl << "Singular value matrix:" << std::endl; + std::cout << S << std::endl << std::endl; + + Eigen::MatrixXd B (NROWS, MCOLS); // N x M + Eigen::MatrixXd V (MCOLS, MCOLS); // M x M + Eigen::MatrixXd VT(MCOLS, MCOLS); // Transpose of V + Eigen::MatrixXd U (NROWS, MCOLS); // N x M + + V=svd.matrixV(); + U=svd.matrixU(); + VT=(svd.matrixV()).transpose(); + + B = U*S*VT; // Recompose the matrix from the decomposition + + std::cout << "Reconstructed data matrix" << std::endl; + std::cout << B << std::endl << std::endl; + + Eigen::MatrixXd C(NROWS,MCOLS); + C = A - B; + std::cout << "Difference with input data" << std::endl; + std::cout << C << std::endl; + std::cout << "Norm: " << C.norm() << std::endl << std::endl; + + Eigen::MatrixXd D(MCOLS, MCOLS); + D = ((svd.matrixU()).transpose())*U; + std::cout << "Left vector orthogonality" << std::endl; + std::cout << D << std::endl << std::endl; + + D = ( (svd.matrixV()).transpose() ) * svd.matrixV(); + std::cout << "Right vector orthogonality" << std::endl; + std::cout << D << std::endl << std::endl; + + // Sign correction test + { + auto U = svd.matrixU(); + auto S = svd.singularValues(); + auto V = svd.matrixV(); + + MSSA::SvdSignChoice(A, U, S, V, true); + + std::cout << "Corrected left singular vectors:" + << std::endl << U << std::endl << std::endl; + + std::cout << "Corrected right singular vectors:" + << std::endl << V << std::endl << std::endl; + + Eigen::MatrixXd C = A - U * S.asDiagonal() * V.transpose(); + std::cout << "Difference with input data" << std::endl; + std::cout << C << std::endl; + std::cout << "Norm: " << C.norm() << std::endl; + + if (C.norm() > 1.0e-10) { + + for (int k=0; k + svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV); + + std::cout << "Its singular values are:" << std::endl + << svd.singularValues() << std::endl << std::endl; + + std::cout << "Its left singular vectors are the columns of the thin U matrix:" + << std::endl << svd.matrixU() << std::endl << std::endl; + + std::cout << "Its right singular vectors are the columns of the thin V matrix:" + << std::endl << svd.matrixV() << std::endl << std::endl; + + // Test that this SVD is doing what we expect. + // + Eigen::MatrixXd S = svd.singularValues().asDiagonal(); + + std::cout << std::endl << "Singular value matrix:" << std::endl; + std::cout << S << std::endl << std::endl; + + Eigen::MatrixXd B (NROWS, MCOLS); // N x M + Eigen::MatrixXd V (MCOLS, MCOLS); // M x M + Eigen::MatrixXd VT(MCOLS, MCOLS); // Transpose of V + Eigen::MatrixXd U (NROWS, MCOLS); // N x M + + V=svd.matrixV(); + U=svd.matrixU(); + VT=(svd.matrixV()).transpose(); + + B = U*S*VT; // Recompose the matrix from the decomposition + + std::cout << "Reconstructed data matrix" << std::endl; + std::cout << B << std::endl << std::endl; + + Eigen::MatrixXd C(NROWS,MCOLS); + C = A - B; + std::cout << "Difference with input data" << std::endl; + std::cout << C << std::endl; + std::cout << "Norm: " << C.norm() << std::endl << std::endl; + + Eigen::MatrixXd D(MCOLS, MCOLS); + D = ((svd.matrixU()).transpose())*U; + std::cout << "Left vector orthogonality" << std::endl; + std::cout << D << std::endl << std::endl; + + D = ( (svd.matrixV()).transpose() ) * svd.matrixV(); + std::cout << "Right vector orthogonality" << std::endl; + std::cout << D << std::endl << std::endl; + + // Sign correction test + { + auto U = svd.matrixU(); + auto S = svd.singularValues(); + auto V = svd.matrixV(); + + MSSA::SvdSignChoice(A, U, S, V, true); + + std::cout << "Corrected left singular vectors:" + << std::endl << U << std::endl << std::endl; + + std::cout << "Corrected right singular vectors:" + << std::endl << V << std::endl << std::endl; + + Eigen::MatrixXd C = A - U * S.asDiagonal() * V.transpose(); + std::cout << "Difference with input data" << std::endl; + std::cout << C << std::endl; + std::cout << "Norm: " << C.norm() << std::endl; + + if (C.norm() > 1.0e-10) { + + for (int k=0; k