Commit 399a3039 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

restructuring

parent 8a1c9185
Pipeline #118324 failed with stages
in 8 minutes and 49 seconds
......@@ -54,7 +54,6 @@ add_subdirectory(mlmc/src)
set(MLMC_LIBRARIES MLMC sprng fftw3 m ${MPP_LIBRARIES})
set(MLMC_TEST_LIBRARIES MLMC sprng fftw3 m ${MPP_TEST_LIBRARIES})
# Executables
add_executable(MLMC-M++ mlmc/src/Main.cpp)
......
include_directories(pdesolver)
add_subdirectory(pdesolver)
include_directories(problems)
add_subdirectory(problems)
include(main/CMakeLists.txt)
include(montecarlo/CMakeLists.txt)
include(problem/CMakeLists.txt)
include(generators/CMakeLists.txt)
add_library(MLMC STATIC ${assemble} ${main}
${montecarlo} ${problem} ${generators})
......
......@@ -6,131 +6,6 @@
using namespace std;
pair<double, double> linear_fit(vector<double> &x, vector<double> &y) {
long xSize = x.size();
double xSum = 0, ySum = 0, xxSum = 0, xySum = 0, slope, intercept;
for (long i = 0; i < xSize; i++) {
xSum += x[i];
ySum += y[i];
xxSum += x[i] * x[i];
xySum += x[i] * y[i];
}
slope = ((double) xSize * xySum - xSum * ySum) /
((double) xSize * xxSum - xSum * xSum);
intercept = (ySum - slope * xSum) / ((double) xSize);
pair<double, double> res;
res.first = slope;
res.second = intercept;
return res;
}
vector<double> linspace(const double &start, const double &end, int num) {
vector<double> linspaced;
if (num == 0) { return linspaced; }
if (num == 1) {
linspaced.push_back(start);
return linspaced;
}
double delta = (end - start) / (num - 1);
for (int i = 0; i < num - 1; i++) {
linspaced.push_back(start + delta * i);
}
linspaced.push_back(end);
return linspaced;
}
void fft(int N, fftw_complex *in, fftw_complex *out) {
fftw_plan plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(plan);
fftw_destroy_plan(plan);
fftw_cleanup();
}
void fft(vector<complex<double>> &in, vector<complex<double>> &out) {
int N = in.size();
fftw_complex in_arr[N], out_arr[N];
for (int i = 0; i < N; i++) {
in_arr[i][REAL] = in[i].real();
in_arr[i][IMAG] = in[i].imag();
}
fft(N, in_arr, out_arr);
out.resize(N);
for (int i = 0; i < N; i++) {
out[i] = {out_arr[i][REAL], out_arr[i][IMAG]};
}
}
void fft(const vector<double> &in, vector<double> &out) {
int N = in.size();
vector<complex<double>> in_complex;
vector<complex<double>> out_complex;
in_complex.resize(N);
out_complex.resize(N);
for (int i = 0; i < N; i++) {
in_complex[i] = {in[i], 0.0};
}
fft(in_complex, out_complex);
out.resize(N);
for (int i = 0; i < N; i++) {
out[i] = out_complex[i].real();
}
}
void fft2(int N2, int N1, fftw_complex *in, fftw_complex *out) {
fftw_plan plan = fftw_plan_dft_2d(N2, N1, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(plan);
fftw_destroy_plan(plan);
fftw_cleanup();
}
void fft2(const vector<vector<complex<double>>> &in,
vector<vector<complex<double>>> &out) {
int N2 = in.size();
int N1 = in[0].size();
auto in_arr = new fftw_complex[N2 * N1];
auto out_arr = new fftw_complex[N2 * N1];
for (int i = 0; i < N2; i++) {
for (int j = 0; j < N1; j++) {
in_arr[i * N2 + j][REAL] = in[i][j].real();
in_arr[i * N2 + j][IMAG] = in[i][j].imag();
}
}
fft2(N2, N1, in_arr, out_arr);
out.resize(N2);
for (int i = 0; i < N2; i++) {
out[i].resize(N1);
for (int j = 0; j < N1; j++) {
out[i][j] = {out_arr[i * N2 + j][REAL], out_arr[i * N2 + j][IMAG]};
}
}
delete[] in_arr;
delete[] out_arr;
}
void fft2(const vector<vector<double>> &in, vector<vector<double>> &out) {
int N2 = in.size();
int N1 = in[0].size();
vector<vector<complex<double>>> in_complex;
vector<vector<complex<double>>> out_complex;
in_complex.resize(N2);
out_complex.resize(N2);
for (int i = 0; i < N2; i++) {
in_complex[i].resize(N1);
out_complex[i].resize(N2);
for (int j = 0; j < N1; j++) {
in_complex[i][j] = {in[i][j], 0.0};
}
}
fft2(in_complex, out_complex);
out.resize(N2);
for (int i = 0; i < N2; i++) {
out[i].resize(N1);
for (int j = 0; j < N1; j++) {
out[i][j] = out_complex[i][j].real();
}
}
}
double calc_mean(const vector<double> &vec) {
double mean_ = 0.0;
for (auto const &value: vec)
......
......@@ -4,46 +4,7 @@
#include <utility>
#include <map>
#include <vector>
#include <complex>
#include <fftw3.h>
std::pair<double, double> linear_fit(std::vector<double> &x, std::vector<double> &y);
std::vector<double> linspace(const double &start,
const double &end, int num);
void fft(int N, fftw_complex *in, fftw_complex *out);
void fft(std::vector<std::complex<double>> &in,
std::vector<std::complex<double>> &out);
void fft(const std::vector<double> &in,
std::vector<double> &out);
void fft2(int N2, int N1, fftw_complex *in, fftw_complex *out);
void fft2(const std::vector<std::vector<std::complex<double>>> &in,
std::vector<std::vector<std::complex<double>>> &out);
void fft2(const std::vector<std::vector<double>> &in,
std::vector<std::vector<double>> &out);
double calc_mean(const std::vector<double> &vec);
double calc_mean(const std::vector<std::vector<double>> &vec);
double calc_mean(const std::vector<int> &vec);
double calc_mean(const std::vector<std::vector<int>> &vec);
double calc_var(const std::vector<double> &vec);
double calc_var(const std::vector<std::vector<double>> &vec);
double calc_var(const std::vector<int> &vec);
double calc_var(const std::vector<std::vector<int>> &vec);
#include <math.h>
template<typename T>
std::string vec2str(std::vector<T> vec);
......
......@@ -24,7 +24,7 @@ double Exponents::ComputeAlpha(const AveragesMap &avgs) {
avgsY.push_back(-log2(_avgsY[i]));
}
return linear_fit(levels, avgsY).first;
return linearFit(levels, avgsY).first;
}
double Exponents::ComputeBeta(const VariancesMap &vars) {
......@@ -37,7 +37,7 @@ double Exponents::ComputeBeta(const VariancesMap &vars) {
varsY.push_back(-log2(_varsY[i]));
}
return linear_fit(levels, varsY).first;
return linearFit(levels, varsY).first;
}
double Exponents::ComputeGamma(const AveragesMap &avgs) {
......@@ -50,5 +50,23 @@ double Exponents::ComputeGamma(const AveragesMap &avgs) {
avgsCost.push_back(log2(_avgsCost[i]));
}
return linear_fit(levels, avgsCost).first;
return linearFit(levels, avgsCost).first;
}
std::pair<double, double> Exponents::linearFit(vector<double> &x, vector<double> &y) {
long xSize = x.size();
double xSum = 0, ySum = 0, xxSum = 0, xySum = 0, slope, intercept;
for (long i = 0; i < xSize; i++) {
xSum += x[i];
ySum += y[i];
xxSum += x[i] * x[i];
xySum += x[i] * y[i];
}
slope = ((double) xSize * xySum - xSum * ySum) /
((double) xSize * xxSum - xSum * xSum);
intercept = (ySum - slope * xSum) / ((double) xSize);
std::pair<double, double> res;
res.first = slope;
res.second = intercept;
return res;
}
......@@ -30,6 +30,9 @@ struct Exponents {
double ComputeBeta(const VariancesMap &vars);
double ComputeGamma(const AveragesMap &avgs);
std::pair<double, double> linearFit(std::vector<double> &x, std::vector<double> &y);
};
#endif //EXPONENTS_HPP
......@@ -8,4 +8,4 @@ add_library(PDESOLVER STATIC
assembling/PGReactionAssemble.cpp
assembling/DGReactionAssemble.cpp
)
target_link_libraries(PDESOLVER ${MPP_LIBRARIES})
\ No newline at end of file
target_link_libraries(PDESOLVER PROBLEMS ${MPP_LIBRARIES})
\ No newline at end of file
......@@ -3,7 +3,7 @@
#include "Assemble.hpp"
#include "discretization/IDiscretization.hpp"
#include "problem/StochasticEllipticProblem.hpp"
#include "problems/StochasticEllipticProblem.hpp"
typedef std::pair<double, double> FluxPair;
......
......@@ -4,7 +4,7 @@
#include"Assemble.hpp"
#include "Plot.hpp"
#include "discretization/IDiscretization.hpp"
#include "problem/StochasticTransportProblem.hpp"
#include "problems/StochasticTransportProblem.hpp"
typedef std::pair<double, double> RatePair;
......
#ifndef _DG_REACTION_H_
#define _DG_REACTION_H_
#include "problem/StochasticReactionProblem.hpp"
#include "problems/StochasticReactionProblem.hpp"
#include "pdesolver/IStochasticReactionAssemble.hpp"
#include "utility/ctools.hpp"
#include "Plot.hpp"
......
......@@ -2,7 +2,7 @@
#define TUTORIAL_DGTRANSPORT_HPP
#include "pdesolver/IStochasticTransportAssemble.hpp"
#include "problem/StochasticTransportProblem.hpp"
#include "problems/StochasticTransportProblem.hpp"
#include "discretization/DGDiscretization.hpp"
#include "TimeSeries.hpp"
#include "utility/ctools.hpp"
......
......@@ -5,7 +5,7 @@
#include "elements/Elements.hpp"
#include "Plot.hpp"
#include "utility/ctools.hpp"
#include "problem/StochasticReactionProblem.hpp"
#include "problems/StochasticReactionProblem.hpp"
#include "pdesolver/IStochasticReactionAssemble.hpp"
......
set(problem
problem/StochasticProblem.cpp
problem/StochasticEllipticProblem.cpp
problem/StochasticTransportProblem.cpp
)
\ No newline at end of file
set(generators
add_library(PROBLEMS STATIC
IStochasticProblem.cpp
StochasticEllipticProblem.cpp
StochasticTransportProblem.cpp
generators/CirculantEmbedding.cpp
generators/HybridFluxGenerator.cpp
generators/RandomNumberManager.cpp
)
\ No newline at end of file
)
target_link_libraries(PROBLEMS fftw3 ${MLMC_LIBRARIES})
......@@ -12,7 +12,6 @@ SampleGenerator *IStochasticProblem::createGenerator(GeneratorName genName,
// if (genName == "NormalRandomNumberGenerator")
// return new NormalRandomNumberGenerator(meshes);
if (genName == "CirculantEmbedding")
return new MultilevelCirculantEmbedding(meshes);
if (genName == "HybridFluxGenerator")
......
#include "problem/StochasticEllipticProblem.hpp"
#include "problems/StochasticEllipticProblem.hpp"
IStochasticEllipticProblem *
......
#ifndef STOCHASTICELLIPTICPROBLEM_HPP
#define STOCHASTICELLIPTICPROBLEM_HPP
#include "problem/IStochasticProblem.hpp"
#include "problems/IStochasticProblem.hpp"
class IStochasticEllipticProblem : public IStochasticProblem {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment