Commit 67913b57 authored by Tim Buchholz's avatar Tim Buchholz
Browse files

rollback 18_02_20

parent cac275a5
......@@ -6,6 +6,3 @@ Todo.txt
cmake-build-debug/
notebooks/__pycache__/
notebooks/.idea/
*/__pycache__/*
*/.ipynb_checkpoints/*
......@@ -2,37 +2,67 @@ cmake_minimum_required(VERSION 3.5.1)
project(MLMC)
#---------------------------------------------------------------------------------------#
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O0")
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O1")
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O2")
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O3")
# Google test (Placed here to not compile it with mpiCC)
add_subdirectory(mpp/googletest)
include_directories(${PROJECT_BINARY_DIR}/mpp/googletest/googletest/include)
include_directories(${PROJECT_BINARY_DIR}/mpp/googlemock/googletest/include)
set(GTEST_LIB gtest gtest_main)
#---------------------------------------------------------------------------------------#
option(NO_DEPRECATED "Suppress deprecated warnings" OFF)
# MPI
find_package(MPI REQUIRED)
include_directories(${MPI_INCLUDE_PATH})
option(BUILD_TESTS "Build test executables" ON)
set(CMAKE_C_COMPILER ${MPI_C_COMPILER})
set(CMAKE_CXX_COMPILER ${MPI_CXX_COMPILER})
option(USE_CXSC "Use interval arithmetic library cxsc" OFF)
option(PROBLEM_NO_TIME "Time independent problem" OFF)
option(PROBLEM_1D "1 dimensional problem" OFF)
option(PROBLEM_2D "2 dimensional problem" OFF)
option(USE_SUPERLU30 "Use SuperLU 3.0 (old version)" OFF)
option(AFFINE_LINEAR_TRAFO "Only affine linear transformations" OFF)
set(CMAKE_CXX_COMPILE_FLAGS ${CMAKE_CXX_COMPILE_FLAGS} ${MPI_COMPILE_FLAGS})
set(CMAKE_CXX_LINK_FLAGS ${CMAKE_CXX_LINK_FLAGS} ${MPI_LINK_FLAGS})
set(CMAKE_BUILD_TYPE distribution)
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-Ofast")
set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O0")
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O1")
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O2")
#set(CMAKE_CXX_FLAGS_DISTRIBUTION "-O3")
set(CMAKE_CXX_FLAGS "-fPIC -g -Wno-deprecated -std=gnu++11")
#set(CMAKE_CXX_FLAGS "-fPIC -g -Wno-deprecated -Wall -std=gnu++11")
#---------------------------------------------------------------------------------------#
# Set path for mpp
set(PROJECT_MPP_DIR ${PROJECT_SOURCE_DIR}/mpp)
# Blas Lapack
#find_package(FFTW3 REQUIRED) (Can't be found but still seems to link in)
find_package(BLAS REQUIRED)
find_package(LAPACK REQUIRED)
#---------------------------------------------------------------------------------------#
# Include settings for mpp
include(${PROJECT_MPP_DIR}/CMakeListsMpp.inc)
# Manage folder structure build folder
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/data/dual)
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/data/gp)
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/data/py)
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/data/vtk)
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/log)
#---------------------------------------------------------------------------------------#
# Link libraries
#find_package(FFTW3 REQUIRED) (Can't be found but still seems to link in)
link_directories(${PROJECT_SOURCE_DIR}/sprng5/lib)
link_directories(${PROJECT_SOURCE_DIR}/mpp/superlu/lib)
if (UNIX AND NOT APPLE)
set(SUPERLU superlu_4.3)
elseif (APPLE)
set(SUPERLU superlu_4.3_mac10.10)
elseif (WIN32)
set(SUPERLU superlu_4.3_cygwin)
endif ()
#---------------------------------------------------------------------------------------#
# Include directories
include_directories(${PROJECT_SOURCE_DIR}/sprng5/include)
include_directories(${PROJECT_SOURCE_DIR}/mlmc/src)
include_directories(${PROJECT_SOURCE_DIR}/mpp/src)
#---------------------------------------------------------------------------------------#
# Subdirectories
add_subdirectory(mpp)
add_subdirectory(mlmc/src)
#---------------------------------------------------------------------------------------#
......@@ -51,4 +81,4 @@ target_link_libraries(MLMC-M++ MLMC sprng SRC LIB_PS ${SUPERLU} blas lapack fftw
#target_link_libraries(TestCirculantEmbedding MLMC sprng SRC LIB_PS ${SUPERLU}
# blas lapack fftw3 m ${GTEST_LIB})
#target_link_libraries(TestRNManager MLMC sprng fftw3 m ${GTEST_LIB})
#---------------------------------------------------------------------------------------#
\ No newline at end of file
#---------------------------------------------------------------------------------------#
## Getting Started
This project combines MLMC methods and the parallel PDE solving software M++. This guide
walks you through the installation process.
walks you through the installation process and presents some example computations.
### Prerequisites
* This version of M++ uses CMake (https://cmake.org/download/) as building tool.
......@@ -49,8 +49,8 @@ Transform in the West http://fftw.org/). To install it run:
```cmake ..```
```make -j```
### MLMC-M++ Experiments
You can find jupyter notebooks with experiments in the
in the ```notebooks``` directory.
### MLMC-M++ Exercises
You can find a PDF with exercises in the ```doc``` directory
and a complementing jupyter notebook in the ```notebook``` directory.
loadconf = mlmc.conf
loadconf = mlmc_transport.conf
loadconf = mlmc.conf
\ No newline at end of file
#Model = LagrangeFEM
#Model = MixedFEM
#Model = HybridFEM
#Model = DGFEM
Model = DGTransport
Overlap = dG1
flux_alpha = 1 #Upwind: 1, Central: 0
#sign = 1
#penalty = 20
#Problem = StochasticLaplace1D
#Problem = StochasticPollution1D
#Problem = DeterministicPollution1D
#Problem = StochasticLaplace2D
#Problem = StochasticPollution2D
#Problem = DeterministicPollution2D
#Problem = StochasticPollutionCosHat2D
#Problem = DeterministicPollutionCosHat2D
#Problem = DeterministicPollutionCosHatSquare500
Problem = StochasticPollutionMollifiedBar2D
StochasticField = LogNormal
#StochasticField = Gaussian
......@@ -10,12 +32,20 @@ initSampleAmount = 12, 6, 3
epsilon = 0.01
mcOnly = false
epsilonList = 0.01, 0.005, 0.001
epsilon_lst = 0.01, 0.005, 0.001, 0.0005
maxLevel = 7
plevel = 3
degree = 2
#functional = L2
#functional = Energy
#funcitonal = H1
#functional = Outflow
#functional = Goal
functional = Mass
MainVerbose = 1
MCVerbose = 2
MLMCVerbose = 2
PDEVerbose = 1
......@@ -25,6 +55,8 @@ GeneratorPlotting = 4
MCPlotting = 3
MLMCPlotting = 2
enablePython = 1
# Stochastic Field
mean = 1.0
sigma = 1.0
......@@ -55,3 +87,20 @@ LinearVerbose = -1
MuteLevel = -1
NewtonVerbose = -1
DebugLevel = -1
#Transport
scaling = 8
TimeSeries = uniform
startTime = 0.0
endTime = 1.0
Kmax = 250
Keps = 1e-5
#rkorder = 0 #Exponential integrator
#rkorder = 1 #Explicit Euler
#rkorder = 2 #Heun
#rkorder = 2 #Runge
rkorder = -2 #Implicit MP-rule
#rkorder = 1000 #Polynomial Krylov-Arnoldi method
gamma = 0.01
\ No newline at end of file
......@@ -3,8 +3,7 @@ set(MLMC_SRC
Utils.C
MLMCMain.C
MultilevelMonteCarlo.C
utils/IndentedLogger.cpp
mc/MonteCarlo.cpp
MonteCarloLogger.C
mc/MonteCarloElliptic.C
mc/MonteCarloTransport.C
assemble/EllipticAssemble.C
......@@ -13,6 +12,6 @@ set(MLMC_SRC
assemble/DGEllipticAssemble.C
stochastics/RNManager.C
stochastics/CirculantEmbedding.C
problem/StochasticEllipticProblem.h stochastics/SampleGenerator.h stochastics/HybridFluxGenerator.C stochastics/HybridFluxGenerator.h problem/StochasticProblemFactory.hpp mc/MonteCarlo.cpp utils/MultilevelPlotter.cpp utils/MultilevelPlotter.hpp utils/MatrixGraphContainer.cpp utils/MatrixGraphContainer.hpp)
problem/StochasticEllipticProblem.h stochastics/SampleGenerator.h stochastics/HybridFluxGenerator.C stochastics/HybridFluxGenerator.h problem/StochasticProblemFactory.hpp)
add_library(MLMC STATIC ${MLMC_SRC})
......@@ -11,10 +11,7 @@ void MLMCMain::initialize() {
problem = StochasticProblemFactory::CreateStochasticProblemFrom(problemName);
assemble = getAssemble();
plots = new MultilevelPlotter(*meshes);
// auto mgContainer = new MatrixGraphContainer(*meshes);
mlmc = new MultilevelMonteCarlo(initLevels, initSampleAmount, meshes, *plots,
mlmc = new MultilevelMonteCarlo(initLevels, initSampleAmount, meshes,
stochasticField, assemble);
}
......@@ -25,7 +22,7 @@ Meshes *MLMCMain::getMeshes() {
return new Meshes("UnitSquare", pLevel, maxLevel);
if (problemName.find("Square500") != string::npos)
return new Meshes("Square500", pLevel, maxLevel);
Exit("\nMesh not found in " + problemName + "\n")
Exit("\nMesh not found\n")
}
Discretization *MLMCMain::getDiscretization() {
......@@ -51,7 +48,7 @@ Discretization *MLMCMain::getDiscretization() {
return nullptr;
if (modelName == "DGTransport")
return new DGDiscretization(new DGDoF(degree));
Exit("\nDiscretization not found in " + modelName + "\n")
Exit("\nDiscretization not implemented yet\n")
}
StochasticField *MLMCMain::getStochasticField() {
......@@ -89,10 +86,9 @@ Assemble *MLMCMain::getAssemble() {
}
void MLMCMain::Run() {
printInfo();
mlmc->PrintInfo();
printParameters();
meshes->PrintInfo();
// assemble->PrintInfo();
if (experimentName == "ConvergenceTest") {
headConvergenceTest();
......@@ -100,6 +96,9 @@ void MLMCMain::Run() {
mlmc->ShowMCTable();
mlmc->ShowKurtosisWarning();
mlmc->ShowExponentResults();
if (enablePython == 1) {
system("python3 ../tools/plot_statistics.py ConvergenceTest");
}
} else if (experimentName == "MLMCExperiment") {
headMLMCExperiment();
mlmc->Method(epsilon);
......@@ -107,24 +106,35 @@ void MLMCMain::Run() {
mlmc->ShowResultsMLMCRun(epsilon);
mlmc->ShowKurtosisWarning();
mlmc->ShowExponentResults();
if (enablePython == 1 && PPM->master()) {
system("python3 ../tools/plot_statistics.py MLMCExperiment");
system("python3 ../tools/plot_mlmc.py MLMCExperiment");
}
} else if (experimentName == "MLMCOverEpsilon") {
headMLMCOverEpsilon();
mlmc->Method(epsilonList);
mlmc->ShowMCTable();
mlmc->ShowResultsOverEpsilon(epsilonList);
if (enablePython == 1 && PPM->master()) {
system("python3 ../tools/plot_statistics.py MLMCOverEpsilon");
system("python3 ../tools/plot_mlmc.py MLMCOverEpsilon");
}
}
}
void MLMCMain::printInfo() {
Vout(1) << endl
<< "Main Info:" << endl
<< " Problem: " << problemName << endl
<< " Model: " << modelName << endl
<< " StochField: " << fieldName << endl
<< " Experiment: " << experimentName << endl;
void MLMCMain::printParameters() {
mout << endl;
mout << "Problem \t=\t" << problemName << endl;
mout << "Model \t=\t" << modelName << endl;
mout << "Experiment\t=\t" << experimentName << endl;
mout << "StochField\t=\t" << fieldName << endl;
mout << "initLevels\t=\t" << vec2str(initLevels) << endl;
mout << "initSample\t=\t" << vec2str(initSampleAmount) << endl;
if (experimentName == "MLMCExperiment")
Vout(1) << " epsilon: " << epsilon << endl;
mout << "eps \t=\t" << epsilon << endl;
if (experimentName == "MLMCOverEpsilon")
Vout(1) << " epsilonList: " << vec2str(epsilonList) << endl;
Vout(1) << endl;
mout << "eps_lst \t=\t" << vec2str(epsilonList) << endl;
mout << "maxLevel \t=\t" << maxLevel << endl;
mout << "pLevel \t=\t" << pLevel << endl;
mout << "enablePy \t=\t" << enablePython << endl;
}
......@@ -6,17 +6,14 @@
#include "MultilevelMonteCarlo.h"
#include "discretization/DGDiscretization.h"
#include "stochastics/StochasticField.h"
#include "utils/MultilevelPlotter.hpp"
#include "utils/MatrixGraphContainer.hpp"
class MLMCMain {
public:
int verbose = 0;
string problemName = "", modelName = "";
string fieldName = "", experimentName = "";
int enablePython = 1;
double epsilon = 0.1;
int pLevel = 0;
int maxLevel = 8;
......@@ -27,7 +24,6 @@ public:
vector<double> epsilonList = {0.1, 0.05, 0.01};
Meshes *meshes = nullptr;
MultilevelPlotter *plots = nullptr;
Discretization *disc = nullptr;
StochasticField *stochasticField = nullptr;
......@@ -36,8 +32,7 @@ public:
MultilevelMonteCarlo *mlmc = nullptr;
MLMCMain() {
config.get("MainVerbose", verbose);
config.get("enablePython", enablePython);
config.get("epsilon", epsilon);
config.get("plevel", pLevel);
config.get("maxLevel", maxLevel);
......@@ -55,13 +50,6 @@ public:
if (maxLevel < initLevels.back()) Exit(
"\nLast element of initLevels has to be smaller or equal to maxLevel.\n")
int uniformSampleAmount;
config.get("uniformSampleAmount", uniformSampleAmount);
if (experimentName == "ConvergenceTest")
fill(initSampleAmount.begin(),
initSampleAmount.end(),
uniformSampleAmount);
initialize();
}
......@@ -89,7 +77,7 @@ private:
Assemble *getAssemble();
void printInfo();
void printParameters();
static void headConvergenceTest() {
mout << endl << "**********************************************************"
......
#include "MonteCarloLogger.h"
IndentedLogger *IndentedLogger::instance = nullptr;
#ifndef MLMC__MULTILEVELMONTECARLO__HPP
#define MLMC__MULTILEVELMONTECARLO__HPP
#include "m++.h"
class IndentedLogger {
static IndentedLogger *instance;
string indent;
string innerIndent;
const string defaultIndent = " ";
vector<Date> timestamps;
IndentedLogger() {
indent = "";
}
public:
static IndentedLogger *GetInstance() {
if (!instance)
instance = new IndentedLogger;
return instance;
}
int DecreaseIndent() {
indent = indent.erase(0, defaultIndent.length());
}
void IncreaseIndent() {
indent = indent.append(defaultIndent);
}
string GetIndent() {
return indent;
}
void StartMethod(const string& msg, int verbose) {
Date start;
timestamps.push_back(start);
vout (1) << indent << "<" << msg << ">" << endl;
IncreaseIndent();
}
void EndMethod(int verbose) {
DecreaseIndent();
vout (1) << indent << "<End after " << Date() - timestamps.back() << ">" << endl;
timestamps.erase(timestamps.end());
}
void InnerMsg(const string& msg, int verbose) {
vout (1) << indent << msg << endl;
}
void LogMsgFlush(const string &msg, int verbose) {
vout(1) << "\r" << indent << msg << flush;
}
};
#endif //MLMC__MULTILEVELMONTECARLO__HPP
......@@ -8,11 +8,11 @@ MonteCarlo *MultilevelMonteCarlo::getMonteCarlo(int l, int dM, bool onlyFineLeve
auto transportAssemble = dynamic_cast<DGTransportAssemble *>(assemble);
if (ellipticAssemble != nullptr)
return new MonteCarloElliptic(l, dM, onlyFineLevel, meshes,
stochasticField, ellipticAssemble, plots);
stochasticField, ellipticAssemble);
if (transportAssemble != nullptr)
return new MonteCarloTransport(l, dM, onlyFineLevel, meshes, plots,
return new MonteCarloTransport(l, dM, onlyFineLevel, meshes,
stochasticField, transportAssemble);
Exit("Wrong Assembly Class")
return nullptr;
}
void MultilevelMonteCarlo::initializeMapMonteCarlo() {
......@@ -66,7 +66,8 @@ void MultilevelMonteCarlo::Method(const double eps) {
}
void MultilevelMonteCarlo::Method(const vector<double> &epsLst) {
logger->StartMethod("Start MLMC Method over epsilon", verbose);
logger->StartMethod("Start MLMC Method", verbose);
// logger->IncreaseIndent();
for (auto &eps : epsLst) {
initializeMapMonteCarlo();
initializeValues();
......@@ -76,28 +77,29 @@ void MultilevelMonteCarlo::Method(const vector<double> &epsLst) {
valueOverEpsilon.push_back(value);
costOverEpsilon.push_back(cost);
}
// logger->DecreaseIndent();
logger->EndMethod(verbose);
}
void MultilevelMonteCarlo::estimateExponents(bool excludeBaseLevel) {
vector<double> levelList, avgsY, varsY, costs;
vector<double> level_lst, avgs_Y, vars_Y, costs;
auto iter = mapMonteCarlo.begin();
if (excludeBaseLevel)
iter++;
for (; iter != mapMonteCarlo.end(); iter++) {
levelList.push_back((double) iter->first);
avgsY.push_back(-log2(iter->second->avgY));
varsY.push_back(-log2(iter->second->varY));
level_lst.push_back((double) iter->first);
avgs_Y.push_back(-log2(iter->second->avgY));
vars_Y.push_back(-log2(iter->second->varY));
costs.push_back(log2(iter->second->avgCost));
}
pair<double, double> res = linear_fit(levelList, avgsY);
alpha = res.first;
res = linear_fit(levelList, varsY);
beta = res.first;
res = linear_fit(levelList, costs);
gamma = res.first;
vector<double> v = linear_fit(level_lst, avgs_Y);
alpha = v.front();
v = linear_fit(level_lst, vars_Y);
beta = v.front();
v = linear_fit(level_lst, costs);
gamma = v.front();
logger->InnerMsg("alpha: " + to_string(alpha) +
" beta: " + to_string(beta) +
......@@ -105,15 +107,15 @@ void MultilevelMonteCarlo::estimateExponents(bool excludeBaseLevel) {
}
void MultilevelMonteCarlo::updateSampleAmount(const double &eps) {
int optimalM;
int M_optimal;
double factor = 0.0;
for (auto &mc : mapMonteCarlo)
factor += sqrt(mc.second->varY * mc.second->avgCost);
for (auto &mc : mapMonteCarlo) {
optimalM = (int) (ceil(2 * pow(eps, -2) * factor *
M_optimal = (int) (ceil(2 * pow(eps, -2) * factor *
sqrt(mc.second->varY / mc.second->avgCost)));
if (optimalM == 1) optimalM++; // Hack
mapMonteCarlo[mc.first]->dM = optimalM - mapMonteCarlo[mc.first]->M;
if (M_optimal == 1) M_optimal++; // Hack
mapMonteCarlo[mc.first]->dM = M_optimal - mapMonteCarlo[mc.first]->M;
}
string msg = "dM:";
......@@ -173,15 +175,6 @@ void MultilevelMonteCarlo::wrapUpResults(double &val,
}
}
void MultilevelMonteCarlo::PrintInfo() {
Vout(1) << endl
<< "MLMC Info:" << endl
<< " initLevels: " << vec2str(initLevels) << endl
<< " initSample: " << vec2str(initSampleAmount) << endl
<< " maxLevel: " << maxLevel << endl
<< endl;
}
void MultilevelMonteCarlo::ShowMCTableHead() {
mout << endl
<< " l M E[Qf-Qc] E[Qf] V[Qf-Qc]"
......
#ifndef MULTILEVELMONTECARLO_H
#define MULTILEVELMONTECARLO_H
#include "utils/IndentedLogger.hpp"
#include "MonteCarloLogger.h"
#include "mc/MonteCarloElliptic.h"
#include "mc/MonteCarloTransport.h"
......@@ -43,7 +43,6 @@ public:
map<int, MonteCarlo *> mapMonteCarlo;
Meshes *meshes;
MultilevelPlotter &plots;
StochasticField *stochasticField;
Assemble *assemble;
......@@ -59,13 +58,11 @@ public:
MultilevelMonteCarlo(vector<int> &initLevels,
vector<int> &initSampleAmount,
Meshes *meshes,
MultilevelPlotter &plots,
StochasticField *stochasticField,
Assemble *assemble) :
initLevels(initLevels),
initSampleAmount(initSampleAmount),
meshes(meshes),
plots(plots),
stochasticField(stochasticField),
assemble(assemble) {
......@@ -92,8 +89,6 @@ public:
void Method(const vector<double> &epsLst);
void PrintInfo();
static void ShowMCTableHead();
static void ShowResultsMLMCRunHead();
......
#include "Utils.h"
#include "m++.h"
#define REAL 0
#define IMAG 1
pair<double, double> linear_fit(vector<double> &x, vector<double> &y) {
vector<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++) {
......@@ -12,12 +13,11 @@ pair<double, double> linear_fit(vector<double> &x, vector<double> &y) {
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;
slope = (xSize * xySum - xSum * ySum) / (xSize * xxSum - xSum * xSum);
intercept = (ySum - slope * xSum) / xSize;
vector<double> res;
res.push_back(slope);
res.push_back(intercept);
return res;
}
......@@ -80,12 +80,11 @@ void fft2(int N2, int N1, fftw_complex *in, fftw_complex *out) {
fftw_cleanup();
}
void fft2(const vector<vector<complex<double>>> &in,
vector<vector<complex<double>>> &out) {
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];
auto out_arr= new fftw_complex[N2 * N1];