mslang.cpp 4.88 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);
            }

kraus's avatar
kraus committed
74 75 76 77 78 79 80
            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;

81 82 83 84
            mslang::QuadTree tree;
            tree.bb_m = mslang::BoundingBox(llc, urc);
            tree.objects_m.insert(tree.objects_m.end(), baseBlocks.begin(), baseBlocks.end());
            tree.buildUp();
85

86 87
            out.open("quadtree.gpl");
            tree.writeGnuplot(out);
88 89 90 91

            out.close();

            out.open("particles.gpl");
92 93
            gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

94
            unsigned int n = 0;
95 96 97
            IpplTimings::TimerRef mainTimer = IpplTimings::getTimer("mainTimer");
            IpplTimings::startTimer(mainTimer);

98
            for (unsigned int i = 0; i < N; ++ i) {
99 100 101 102
                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);

103 104 105

                if (method == 0) {
                    if (tree.isInside(X)) {
106
                        ++ n;
107 108 109
                        out << std::setw(14) << X[0]
                            << std::setw(14) << X[1]
                            << std::endl;
110 111 112 113 114 115 116 117 118 119
                    }
                } 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;
                        }
120 121 122
                    }
                }
            }
123 124
            IpplTimings::stopTimer(mainTimer);

125
            std::cout << (double)n / N * 100 << " % of particles passed" << std::endl;
126 127 128 129 130 131 132 133

            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);
134 135
            gsl_rng_free(rng);

kraus's avatar
kraus committed
136

137
            for (mslang::Base* func: baseBlocks) {
138 139
                delete func;
            }
kraus's avatar
kraus committed
140
        }
141 142

        out.close();
143
    }
kraus's avatar
kraus committed
144 145 146

    delete fun;

147 148
    return 0;
}