mslang.cpp 5.14 KB
Newer Older
1
#include "Algorithms/Vektor.h"
2
#include "Utilities/MSLang.h"
kraus's avatar
kraus committed
3

4 5
#include "Ippl.h"
#include "Utility/IpplTimings.h"
6

7 8
#include <gsl/gsl_rng.h>

kraus's avatar
kraus committed
9 10
#include <boost/regex.hpp>

11 12
#include <iostream>
#include <string>
13 14
#include <fstream>
#include <streambuf>
15
#include <cstdlib>
16
#include <cmath>
17
#include <set>
18

19
Ippl *ippl;
20
int main(int argc, char *argv[])
21
{
22 23 24 25
    ippl = new Ippl(argc, argv);

    if (argc < 4) {
        std::cout << "please provide the name of the file that contains your code, the method and the number of particles" << std::endl;
26 27
        return 1;
    }
28
    mslang::Function *fun;
29 30 31 32 33

    std::ifstream t(argv[1]);
    std::string str((std::istreambuf_iterator<char>(t)),
                    std::istreambuf_iterator<char>());

34 35 36 37
    const unsigned int method = (atoi(argv[2]) == 0 ? 0: 1);
    std::cout << method << std::endl;
    const unsigned int N = atoi(argv[3]);

38
    // std::string str("repeat( translate(union(rectangle(0.1, 0.1), ellipse(0.1, 0.1)), -0.01, -0.02), 2, 0.1, 0.2)");
39

kraus's avatar
kraus committed
40 41 42 43
    str = boost::regex_replace(str, boost::regex("//.*?\\n"), std::string(""), boost::match_default | boost::format_all);
    str = boost::regex_replace(str, boost::regex("\\s"), std::string(""), boost::match_default | boost::format_all);

     if (parse(str, fun)) {
44
        fun->print(0);
kraus's avatar
kraus committed
45 46
        std::cout << "\n" << std::endl;

47
        std::vector<mslang::Base*> baseBlocks;
kraus's avatar
kraus committed
48 49
        fun->apply(baseBlocks);

50
        std::ofstream out("test.gpl");
51 52 53
        for (mslang::Base* bfun: baseBlocks) {
            // bfun->print(0);
            // std::cout << std::endl;
54
            bfun->computeBoundingBox();
55
            bfun->writeGnuplot(out);
kraus's avatar
kraus committed
56
        }
57
        out.close();
58 59 60

        if (baseBlocks.size() > 0) {
            Vector_t llc, urc;
61 62
            mslang::Base* first = baseBlocks.front();
            const mslang::BoundingBox &bb = first->bb_m;
63 64 65 66 67 68 69
            llc = Vector_t(bb.center_m[0] - 0.5 * bb.width_m,
                           bb.center_m[1] - 0.5 * bb.height_m,
                           0.0);
            urc = Vector_t(bb.center_m[0] + 0.5 * bb.width_m,
                           bb.center_m[1] + 0.5 * bb.height_m,
                           0.0);

70

71
            for (unsigned int i = 1; i < baseBlocks.size(); ++ i) {
72
                const mslang::BoundingBox &bb = baseBlocks[i]->bb_m;
73 74 75 76 77 78
                llc[0] = std::min(llc[0], bb.center_m[0] - 0.5 * bb.width_m);
                llc[1] = std::min(llc[1], bb.center_m[1] - 0.5 * bb.height_m);
                urc[0] = std::max(urc[0], bb.center_m[0] + 0.5 * bb.width_m);
                urc[1] = std::max(urc[1], bb.center_m[1] + 0.5 * bb.height_m);
            }

kraus's avatar
kraus committed
79 80 81 82 83 84 85
            double width = urc[0] - llc[0];
            double height = urc[1] - llc[1];
            llc[0] -= 1e-3 * width;
            urc[0] += 1e-3 * width;
            llc[1] -= 1e-3 * height;
            urc[1] += 1e-3 * width;

86 87 88 89
            mslang::QuadTree tree;
            tree.bb_m = mslang::BoundingBox(llc, urc);
            tree.objects_m.insert(tree.objects_m.end(), baseBlocks.begin(), baseBlocks.end());
            tree.buildUp();
90

91 92
            out.open("quadtree.gpl");
            tree.writeGnuplot(out);
93 94 95 96

            out.close();

            out.open("particles.gpl");
97 98
            gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

99
            unsigned int n = 0;
100 101 102
            IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("mainTimer");
            IpplTimings::startTimer(mainTimer);

103
            for (unsigned int i = 0; i < N; ++ i) {
104 105 106 107
                Vector_t X(0.0);
                X[0] = llc[0] + (urc[0] - llc[0]) * gsl_rng_uniform(rng);
                X[1] = llc[1] + (urc[1] - llc[1]) * gsl_rng_uniform(rng);

108 109 110

                if (method == 0) {
                    if (tree.isInside(X)) {
111
                        ++ n;
112 113 114
                        out << std::setw(14) << X[0]
                            << std::setw(14) << X[1]
                            << std::endl;
115 116 117 118 119 120 121 122 123 124
                    }
                } else {
                    for (mslang::Base* func: baseBlocks) {
                        if (func->isInside(X)) {
                            ++ n;
                            out << std::setw(14) << X[0]
                                << std::setw(14) << X[1]
                                << std::endl;
                            break;
                        }
125 126 127
                    }
                }
            }
128 129
            IpplTimings::stopTimer(mainTimer);

130
            std::cout << (double)n / N * 100 << " % of particles passed" << std::endl;
131 132 133 134 135 136 137 138

            std::map<std::string, unsigned int> characteristicValues;
            characteristicValues.insert(std::make_pair("method", method));
            characteristicValues.insert(std::make_pair("num particles", N));
            characteristicValues.insert(std::make_pair("num base functions", baseBlocks.size()));
            std::stringstream ss;
            ss << "timing__m_=_" << method << "__np_=_" << N << "__nbf_=_" << baseBlocks.size() << ".dat";
            IpplTimings::print(ss.str(), characteristicValues);
139 140
            gsl_rng_free(rng);

kraus's avatar
kraus committed
141

142
            for (mslang::Base* func: baseBlocks) {
143 144
                delete func;
            }
kraus's avatar
kraus committed
145
        }
146 147

        out.close();
148
    }
kraus's avatar
kraus committed
149 150 151

    delete fun;

152 153
    return 0;
}