mslang.cpp 4.65 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>

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

17
Ippl *ippl;
18
int main(int argc, char *argv[])
19
{
20 21 22 23
    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;
24 25
        return 1;
    }
26
    mslang::Function *fun;
27 28 29 30 31

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

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

36
    // 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)");
37 38 39

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

42
        std::vector<mslang::Base*> baseBlocks;
kraus's avatar
kraus committed
43 44
        fun->apply(baseBlocks);

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

        if (baseBlocks.size() > 0) {
            Vector_t llc, urc;
56 57
            mslang::Base* first = baseBlocks.front();
            const mslang::BoundingBox &bb = first->bb_m;
58 59 60 61 62 63 64
            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);

65

66
            for (unsigned int i = 1; i < baseBlocks.size(); ++ i) {
67
                const mslang::BoundingBox &bb = baseBlocks[i]->bb_m;
68 69 70 71 72 73
                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);
            }

74 75 76 77
            mslang::QuadTree tree;
            tree.bb_m = mslang::BoundingBox(llc, urc);
            tree.objects_m.insert(tree.objects_m.end(), baseBlocks.begin(), baseBlocks.end());
            tree.buildUp();
78

79 80
            out.open("quadtree.gpl");
            tree.writeGnuplot(out);
81 82 83 84

            out.close();

            out.open("particles.gpl");
85 86
            gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

87
            unsigned int n = 0;
88 89 90
            IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("mainTimer");
            IpplTimings::startTimer(mainTimer);

91
            for (unsigned int i = 0; i < N; ++ i) {
92 93 94 95
                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);

96 97 98

                if (method == 0) {
                    if (tree.isInside(X)) {
99
                        ++ n;
100 101 102
                        out << std::setw(14) << X[0]
                            << std::setw(14) << X[1]
                            << std::endl;
103 104 105 106 107 108 109 110 111 112
                    }
                } 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;
                        }
113 114 115
                    }
                }
            }
116 117
            IpplTimings::stopTimer(mainTimer);

118
            std::cout << (double)n / N * 100 << " % of particles passed" << std::endl;
119 120 121 122 123 124 125 126

            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);
127 128
            gsl_rng_free(rng);

kraus's avatar
kraus committed
129

130
            for (mslang::Base* func: baseBlocks) {
131 132
                delete func;
            }
kraus's avatar
kraus committed
133
        }
134 135

        out.close();
136
    }
kraus's avatar
kraus committed
137 138 139

    delete fun;

140 141
    return 0;
}