Skip to content
Snippets Groups Projects
Commit a6a2946b authored by Florian Wolf's avatar Florian Wolf
Browse files

initial import from pre-existing code

parents
Branches
No related tags found
No related merge requests found
*~
DrawArea
DrawArea-mpi
#!/bin/bash
#$ -S /bin/bash
#$ -l h_rt=72:00:00
#$ -l h_vmem=6G
#$ -pe openmpi-orte-20 20-1000
#$ -o /work/$USER/$JOB_NAME-$JOB_ID.log
#$ -j y
module load openmpi/gcc
export PSM_RANKS_PER_CONTEXT=4
mpirun -np $NSLOTS \
~/DrawArea/DrawArea-mpi \
-o /work/$USER/$JOB_NAME-$JOB_ID.pgm \
-s 60x18+758+820 -a 100 -m 72 -l 28800 \
~/DrawArea/spain_time_100m.pgm
#!/bin/bash
#$ -S /bin/bash
#$ -l h_rt=72:00:00
#$ -l h_vmem=6G
#$ -binding linear:1
#$ -o /work/$USER/$JOB_NAME-$JOB_ID.log
#$ -j y
~/DrawArea/DrawArea \
-o /work/$USER/$JOB_NAME-$JOB_ID.pgm \
-s 60x18+758+820 -a 100 -m 72 -l 28800 \
~/DrawArea/spain_time_100m.pgm
<?xml version="1.0" encoding="UTF-8" standalone="yes" ?>
<CodeBlocks_project_file>
<FileVersion major="1" minor="6" />
<Project>
<Option title="DrawArea" />
<Option makefile="makefile" />
<Option pch_mode="2" />
<Option compiler="gcc" />
<Build>
<Target title="Release">
<Option output="DrawArea" prefix_auto="1" extension_auto="1" />
<Option object_output="obj/" />
<Option type="1" />
<Option compiler="gcc" />
<Option parameters="-h" />
<Compiler>
<Add option="-O2" />
</Compiler>
<Linker>
<Add option="-s" />
</Linker>
</Target>
</Build>
<Compiler>
<Add option="-O3" />
<Add option="-Wall" />
<Add option="-fexceptions" />
</Compiler>
<Unit filename="src/ClosedList.cpp" />
<Unit filename="src/DrawArea.cpp" />
<Unit filename="src/DrawArea.h" />
<Unit filename="src/Error.h" />
<Unit filename="src/Graph.cpp" />
<Unit filename="src/OpenList.cpp" />
<Unit filename="src/ParseArgs.h" />
<Unit filename="src/WayPoint.cpp" />
<Unit filename="src/pnm.h" />
<Extensions>
<code_completion />
<envvars />
<debugger />
<lib_finder disable_auto="1" />
</Extensions>
</Project>
</CodeBlocks_project_file>
#!/bin/bash
#-----------------------------------------------
#-----------------------------------------------
# job submission loop
#-----------------------------------------------
for i in `/home/wolffl/DrawArea/DrawArea -c 23 -r 12 -l 28800 -m 72 timewalk_100.pgm | sed "s/scope [0-9]* is //"`; do
qsub /home/wolffl/DrawArea/DrawArea_wolffl.sh $i
done
#!/bin/bash
#$ -S /bin/bash
#$ -l h_rt=22:00:00
#$ -l h_vmem=1G
#$ -binding linear:1
#$ -o /work/$USER/$JOB_NAME-$JOB_ID.log
#$ -j y
AREA=$1
/home/wolffl/DrawArea/DrawArea \
-o /work/$USER/area_of_spain_timewalk_100.pgm \
-s $AREA -a 100 -m 72 -l 28800 \
/home/wolffl/DrawArea/timewalk_100.pgm
This diff is collapsed.
makefile 0 → 100644
CC = g++
CFLAGS = -Wall -O3 -c
LFLAGS = -s -lm
PROG = DrawArea
HDR_FILES = DrawArea.h pnm.h ParseArgs.h Error.h
OBJ_FILES = Graph.o WayPoint.o ClosedList.o OpenList.o DrawArea.o
SRC_DIR = ./src
OBJ_DIR = ./obj
HDR = $(HDR_FILES:%=$(SRC_DIR)/%)
OBJ = $(OBJ_FILES:%=$(OBJ_DIR)/%)
$(PROG): $(OBJ) $(HDR)
$(CC) $(LFLAGS) -o $(PROG) $(OBJ)
$(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp object
$(CC) $(CFLAGS) $< -o $@
object:
mkdir -p $(OBJ_DIR)
clean:
rm -rf $(OBJ_DIR) $(PROG)
CC = mpiCC -DDRAWAREA_MPI
CFLAGS = -Wall -O3 -c
LFLAGS = -s -lm
PROG = DrawArea-mpi
HDR_FILES = DrawArea.h pnm.h ParseArgs.h Error.h
OBJ_FILES = Graph.o WayPoint.o ClosedList.o OpenList.o DrawArea.o
SRC_DIR = ./src
OBJ_DIR = ./obj
HDR = $(HDR_FILES:%=$(SRC_DIR)/%)
OBJ = $(OBJ_FILES:%=$(OBJ_DIR)/%)
$(PROG): $(OBJ) $(HDR)
$(CC) $(LDFLAGS) $(LFLAGS) -o $(PROG) $(OBJ)
$(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp object
$(CC) $(CPPFLAGS) $(CFLAGS) $< -o $@
object:
mkdir -p $(OBJ_DIR)
clean:
rm -rf $(OBJ_DIR) $(PROG)
/*********************************************************************
Copyright 2016, Michael Steinert
This file is part of DrawArea (Dijkstra Routing Algorithm Which
Ascertains REachable Areas).
DrawArea is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
DrawArea is distributed WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DrawArea. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#include "DrawArea.h"
// constructor (unneccessary but just for completeness)
ClosedList::ClosedList() {}
// destructor (see constructor)
ClosedList::~ClosedList() {}
// add MapPoint to list
void ClosedList::Push(Node* node)
{
v.insert(node);
}
// get size of list
size_t ClosedList::Size()
{
return v.size();
}
// check if list is empty
bool ClosedList::Empty()
{
return v.empty();
}
// check if MapPoint is contained in list
bool ClosedList::Contains(Node* node)
{
return v.find(node)!=v.end();
}
/*********************************************************************
Copyright 2016, Michael Steinert
This file is part of DrawArea (Dijkstra Routing Algorithm Which
Ascertains REachable Areas).
DrawArea is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
DrawArea is distributed WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DrawArea. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
/*************************************************
TODO-List:
- implement multi-threading
- change temporal distance calculation to
topographic maps (incl. heuristic splitting)
**************************************************/
#ifdef DRAWAREA_MPI
#include "mpi.h"
#define MPI_MSG(msg) "Proc %d " msg, ID
#define MPI_ROOT_PROC 0
#define MPI_STDOUT stderr
#else
#define MPI_MSG(msg) msg
#define MPI_STDOUT stdout
#endif
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include "DrawArea.h"
#include "ParseArgs.h"
#include "Error.h"
#include "pnm.h"
using namespace std;
int ID;
Node* ReadScopeFromMatrix(const char* fileName, Scope& scope);
int SplitHeuristic(const char* inputName, Scope inputScope, Scope* scopeList, size_t nColumns, size_t nRows, float limitTime);
int CalculateAreaMap(const char* inputName, const char* outputName, Scope& scope, float limitTime, float minTime, float areaScaling, int nBinning);
Scope& FitBinning(Scope& scope, int nBinning, int xOffset, int yOffset);
size_t* SplitToEqualParts(size_t nValues, double* values, size_t nParts, size_t *parts=NULL);
double PixelWeight(double time, double limit)
{
return pow(max(1.0-time/limit,0.0),4);
}
char helpText[]="Usage: %s [OPTIONS] -l TIME -m TIME input\n"
"\nMandatories:\n"
" -l, --limit TIME\n"
"\ttravel time limit\n"
" -m, --min-time TIME\n"
"\tminimum time to pass single pixel\n"
"\nOptions:\n"
" -h, --help\n"
"\tprint this message\n"
" -a, --area AREA\n"
"\tarea scaling for output data in pixels / AREA\n"
" -b, --binning N\n"
"\tcalculate only every Nth point in each row and column\n"
" -o, --output FILE\n"
"\tforce output to be written to FILE\n"
" -s, --scope WIDTHxHEIGHT+LEFT+TOP\n"
"\tset scope of interest for map calculation\n"
#ifndef DRAWAREA_MPI
" -c, --cols N\n"
"\tforce map splitting into [N] columns\n"
" -r, --rows N\n"
"\tforce map splitting into N rows\n"
#endif
#ifdef DRAWAREA_THREADS
" -n, --no-calc\n"
"\tdon't perform calculation,\n"
"\tprint only split scopes and prepare output file\n"
" -t, --threads N\n"
"\tsplit map calculation into N threads\n"
"\tIf the number of columns or rows and threads is given, the requested number of threads is ignored.\n"
#endif
"\nTIME and AREA must be positive numbers; N must be a positive integer; WIDTH, HEIGHT, LEFT, and TOP must be integers.\n";
int main(int argc, char* argv[])
{
int error;
int xOffsetBinning;
int yOffsetBinning;
Scope* scopeList;
CmdArgStruct params;
// set error and help message parameters
InitErrorHandling(ExtractAppName(argv[0]),helpText,true);
#ifdef DRAWAREA_MPI
int NoP;
unsigned long yOffset, height;
unsigned long xOffset, width;
unsigned long *heightList, *yOffsetList;
unsigned long *widthList, *xOffsetList;
// init MPI, get process id and total number of processes
if (MPI_Init(&argc,&argv)!=MPI_SUCCESS)
return Error(-1, "MPI_Init failed.\n");
MPI_Comm_size(MPI_COMM_WORLD, &NoP);
MPI_Comm_rank(MPI_COMM_WORLD, &ID);
if (ID==MPI_ROOT_PROC)
{
#endif // DRAWAREA_MPI
// help requested?
if (getCmdOption(argv+1,argv+argc,"-h") || getCmdOption(argv+1,argv+argc,"--help"))
return Help();
// evaluate cmd options
if ((error=params.Evaluate(argc,argv))!=ERROR_SUCCESS)
return error;
// check if input file readable in PNM format
PortableAnyMap inputMap(params.inputName,PNM_MODE_READ);
if (!inputMap.IsOK())
{
return Error(inputMap.GetLastError(),"Input data file %s could not be read!\n%s\n",params.inputName,inputMap.GetLastErrorMessage());
}
params.scope.Fit(inputMap.GetSize());
xOffsetBinning=((inputMap.GetWidth()-1)%params.nBinning)/2;
yOffsetBinning=((inputMap.GetHeight()-1)%params.nBinning)/2;
// check if output file is ready for updating (opens in update mode, has the write type and min size)
// else try to create new file and open for update
Scope outputScope((inputMap.GetWidth()-1)/params.nBinning+1,(inputMap.GetHeight()-1)/params.nBinning+1);
PortableAnyMap outputMap(params.outputName,PNM_MODE_UPDATE);
if (!outputMap.IsOK() || outputMap.GetType()!=PNM_TYPE_P5 || outputMap.GetWidth()<outputScope.Width() || outputMap.GetHeight()<outputScope.Height())
{
Error(ERROR_WARNING,"Output file %s not ready for update. Trying to recreate...\n",params.outputName);
outputMap.Close();
outputMap.Create(params.outputName,outputScope,PNM_TYPE_P5,65535);
outputMap.Open(params.outputName,PNM_MODE_UPDATE);
if (!outputMap.IsOK())
return Error(outputMap.GetLastError(),"Failed to open output file %s.\n%s\n",params.outputName,outputMap.GetLastErrorMessage());
}
inputMap.Close();
outputMap.Close();
#ifdef DRAWAREA_MPI
params.nRows=NoP;
params.nColumns=1;
params.nThreads=NoP;
#endif // DRAWAREA_MPI
// setup height and offset list for heuristic splitting
scopeList=new Scope[params.nThreads];
if (!scopeList)
return Error(-1, MPI_MSG("Unable to allocate memory for heuristic splitting.\n"));
// split map into params.nRows*params.nColumns scopes
// heuristic will split in scopes with approx. equal calculation time
if ((error=SplitHeuristic(params.inputName, params.scope, scopeList, params.nColumns, params.nRows, params.limitTime))!=ERROR_SUCCESS)
return error;
for (int i=0; i<params.nThreads; i++)
FitBinning(scopeList[i],params.nBinning,xOffsetBinning,yOffsetBinning);
#ifdef DRAWAREA_MPI
// prepare scopes for MPI scattering
widthList=new unsigned long[params.nThreads];
xOffsetList=new unsigned long[params.nThreads];
heightList=new unsigned long[params.nThreads];
yOffsetList=new unsigned long[params.nThreads];
if (!widthList || !xOffsetList || !heightList || !yOffsetList)
return Error(-1, MPI_MSG("Unable to allocate memory for mpi scattering of scopes.\n"));
for (int i=0; i<params.nThreads; i++)
{
widthList[i]=scopeList[i].Width();
xOffsetList[i]=scopeList[i].First(COLUMN);
heightList[i]=scopeList[i].Height();
yOffsetList[i]=scopeList[i].First(ROW);
}
}
// send evaluated cmd args to all processes
MPI_Bcast(&params.inputName, 256, MPI_UNSIGNED_CHAR, MPI_ROOT_PROC, MPI_COMM_WORLD);
MPI_Bcast(&params.outputName, 256, MPI_UNSIGNED_CHAR, MPI_ROOT_PROC, MPI_COMM_WORLD);
MPI_Bcast(&params.areaScaling, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD);
MPI_Bcast(&params.limitTime, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD);
MPI_Bcast(&params.minTime, 1, MPI_FLOAT, MPI_ROOT_PROC, MPI_COMM_WORLD);
MPI_Bcast(&params.nBinning, 1, MPI_INT, MPI_ROOT_PROC, MPI_COMM_WORLD);
MPI_Bcast(&xOffsetBinning, 1, MPI_INT, MPI_ROOT_PROC, MPI_COMM_WORLD);
MPI_Bcast(&yOffsetBinning, 1, MPI_INT, MPI_ROOT_PROC, MPI_COMM_WORLD);
// send calculated scopes to all processes
MPI_Scatter(yOffsetList, 1, MPI_UNSIGNED_LONG, &yOffset, 1, MPI_UNSIGNED_LONG, MPI_ROOT_PROC, MPI_COMM_WORLD);
MPI_Scatter(heightList, 1, MPI_UNSIGNED_LONG, &height, 1, MPI_UNSIGNED_LONG, MPI_ROOT_PROC, MPI_COMM_WORLD);
MPI_Scatter(xOffsetList, 1, MPI_UNSIGNED_LONG, &xOffset, 1, MPI_UNSIGNED_LONG, MPI_ROOT_PROC, MPI_COMM_WORLD);
MPI_Scatter(widthList, 1, MPI_UNSIGNED_LONG, &width, 1, MPI_UNSIGNED_LONG, MPI_ROOT_PROC, MPI_COMM_WORLD);
// set process scope
params.scope=Scope(width,height,xOffset,yOffset);
#else // ifdef DRAWAREA_MPI
if (params.nThreads>1)
{
for (int i=0; i<params.nThreads; i++)
printf("scope %d is %ux%u+%u+%u\n",i,scopeList[i].Width(),scopeList[i].Height(),scopeList[i].First(COLUMN),scopeList[i].First(ROW));
#ifdef DRAWAREA_THREADS
if (params.noCalc)
#endif
return 0;
}
#endif // DRAWAREA_MPI
// perform local calculation
clock_t start=clock();
CalculateAreaMap(params.inputName, params.outputName, params.scope, params.limitTime, params.minTime, params.areaScaling, params.nBinning);
fprintf(MPI_STDOUT, MPI_MSG("map calculation took %.2f s\n"),float(clock()-start)/CLOCKS_PER_SEC);
#ifdef DRAWAREA_MPI
MPI_Finalize();
#endif // DRAWAREA_MPI
return 0;
}
// shrink scope to matches binning offsets
Scope& FitBinning(Scope& scope, int nBinning, int xOffset, int yOffset)
{
size_t last;
last=scope.Last(COLUMN)-(nBinning+scope.Last(COLUMN)%nBinning-xOffset)%nBinning;
scope.First(COLUMN)+=(nBinning+xOffset-scope.First(COLUMN)%nBinning)%nBinning;
scope.Width()=(last<scope.First(COLUMN)?0:last-scope.First(COLUMN)+1);
last=scope.Last(ROW)-(nBinning+scope.Last(ROW)%nBinning-yOffset)%nBinning;
scope.First(ROW)+=(nBinning+yOffset-scope.First(ROW)%nBinning)%nBinning;
scope.Height()=(last<scope.First(ROW)?0:last-scope.First(ROW)+1);
return scope;
}
// split a list of nValues values into nParts with approx. equal sum(values)
size_t* SplitToEqualParts(size_t nValues, double* values, size_t nParts, size_t *parts)
{
// get total weight
double totalWeight=0.0;
for (size_t i=0; i<nValues; i++)
totalWeight+=values[i];
// allocate memory for part split
if (!parts)
parts=new size_t[nParts+1];
if (!parts)
return NULL;
size_t part=0;
double partWeightActual=0.0;
double partWeightSetPoint=totalWeight/(nParts-part);
parts[part]=0;
for (size_t i=0; i<nValues && part<nParts-1; i++)
{
// if not last remaining part, part contains at least 1 value,
// and actual part weight plus half value exceeds set point
// switch to new part
if ((part<nParts-1) && (parts[part]<i) && (partWeightActual+values[i]/2>partWeightSetPoint))
{
parts[++part]=i;
totalWeight-=partWeightActual;
partWeightActual=0.0;
partWeightSetPoint=totalWeight/(nParts-part);
}
// add value to part
partWeightActual+=values[i];
}
// set remaining parts to end of value list
while (part<nParts)
parts[++part]=nValues;
return parts;
}
// read input file and split into nColumns x nRows equally weighted scopes
// - the weight of each scope is the sum of the evaluated pixel times to the power of 4
// - the evaluated pixel time is MAX(1.0 - PIXEL_TIME / TIME_LIMIT, 0.0)
int SplitHeuristic(const char* inputName, Scope inputScope, Scope* scopeList, size_t nColumns, size_t nRows, float limitTime)
{
// check if scopeList is initialized
if (!scopeList)
return Error(-1, MPI_MSG("No memory allocated for scope list.\n"));
// open input map
PortableAnyMap inputMap(inputName,PNM_MODE_READ);
if (!inputMap.IsOK())
return Error(inputMap.GetLastError(),"Input data file %s could not be read!\n%s\n", inputName,inputMap.GetLastErrorMessage());
inputScope.Fit(inputMap.GetSize());
// allocate memory for map evaluation
size_t *pixs=new size_t[inputScope.Width()];
double *weights=new double[std::max(inputScope.Height(),inputScope.Width())];
size_t *parts=new size_t[std::max(nRows,nColumns)+1];
if (!pixs || !weights || !parts)
return Error(-1, MPI_MSG("Unable to allocate memory for heuristics.\n"));
// row splitting first:
if (nRows>1)
{
// get weight of each line in map and total map weight
Scope readScope=inputScope;
readScope.Height()=1;
for (size_t i=0; i<inputScope.Height(); i++)
{
weights[i]=0.0;
size_t num=inputMap.ReadScope(readScope,pixs);
for (size_t j=0; j<num; j++)
weights[i]+=PixelWeight(pixs[j],limitTime);
readScope.First(ROW)++;
}
}
SplitToEqualParts(inputScope.Height(),weights,nRows,parts);
for (size_t row=0; row<nRows; row++)
for (size_t col=0; col<nColumns; col++)
scopeList[row*nColumns+col].GetRange(ROW)=Range(parts[row+1]-parts[row],inputScope.First(ROW)+parts[row]);
// loop through all row bundles to split into columns
for (size_t row=0; row<nRows; row++)
{
if (nColumns>1)
{
// get weight of each column in row bundle and total bundle weight
Scope readScope=inputScope;
readScope.Height()=1;
readScope.First(ROW)=scopeList[row*nColumns].First(ROW);
for (size_t i=0; i<inputScope.Width(); i++)
weights[i]=0.0;
for (size_t i=0; i<scopeList[row*nColumns].Height(); i++)
{
size_t num=inputMap.ReadScope(readScope,pixs);
for (size_t j=0; j<num; j++)
weights[j]+=PixelWeight(pixs[j],limitTime);
readScope.First(ROW)++;
}
}
// get splitting and save in scopes
SplitToEqualParts(inputScope.Width(),weights,nColumns,parts);
for (size_t col=0; col<nColumns; col++)
scopeList[row*nColumns+col].GetRange(COLUMN)=Range(parts[col+1]-parts[col],parts[col]+inputScope.First(COLUMN));
}
inputMap.Close();
delete [] pixs;
delete [] weights;
delete [] parts;
return ERROR_SUCCESS;
}
// read Node array scope from map file
// if scope is out of total map boundaries, adjust scope
Node* ReadScopeBinaryMatrix(const char* fileName, Scope& scope)
{
// open input map
PortableAnyMap input(fileName,PNM_MODE_READ);
if (!input.IsOK())
{
Error(input.GetLastError(), MPI_MSG("Unable to open matrix file %s to read data.\n%s\n"), fileName, input.GetLastErrorMessage());
return NULL;
}
// adjust scope to fit into map dimensions
scope.Fit(input.GetSize());
// allocate memory for Node array
Node* nodes=new Node[scope.Size()];
if (!nodes)
{
Error(-1, MPI_MSG("Unable to allocate memory for nodes.\n"));
return NULL;
}
fprintf(MPI_STDOUT, MPI_MSG("reading nodes from scope %u x %u + %u + %u\n"),scope.Width(),scope.Height(),scope.First(COLUMN),scope.First(ROW));
// read nodes scope from input map
input.ReadScope(scope,nodes),scope.Size();
// close input file
input.Close();
return nodes;
}
// calculate area in scope:
// - create maximal reachable node map from
// input map file
// - calculate area for each node in scope
// and write it to output map file
int CalculateAreaMap(const char* inputName, const char* outputName, Scope& scope, float limitTime, float minTime, float areaScaling, int nBinning)
{
// open map file to update pixel data
PortableAnyMap output(outputName,PNM_MODE_UPDATE);
if (!output.IsOK())
return Error(output.GetLastError(), MPI_MSG("Unable to open matrix file %s for data output.\n%s\n"), outputName, output.GetLastErrorMessage());
// adjust scope, if necessary
scope.Fit(output.GetSize()*nBinning);
// enlarge scope by the maximal travel distance in TIME_LIMIT
// to get read scope (all nodes that may be needed)
Scope readScope=scope+int(limitTime/minTime);
// read nodes from binary file
Node* nodes=ReadScopeBinaryMatrix(inputName,readScope);
if (!nodes)
return -1;
// connect nodes
Edge::minTime=minTime;
CreateGraphEdges(nodes,readScope.Width(),readScope.Height());
// adjust scope if necessary
scope.Fit(readScope);
Scope writeScope((scope.Width()-1)/nBinning+1,1,scope.First(COLUMN)/nBinning,scope.First(ROW)/nBinning);
// calculate offset between read nodes and scope nodes
size_t nodesOffset=(scope.First(ROW)-readScope.First(ROW))*readScope.Width()
+ scope.First(COLUMN)-readScope.First(COLUMN);
// allocate memory to write one row at once
size_t* area=new size_t[writeScope.Width()];
if (!area)
return Error(-1, MPI_MSG("No memory allocated for row data.\n"));
fprintf(MPI_STDOUT, MPI_MSG("calculating area in scope %u x %u + %u + %u\n"),scope.Width(),scope.Height(),scope.First(COLUMN),scope.First(ROW));
clock_t start;
// loop all nodes in scope and
// get area that could be reached within TIME_LIMIT
for (size_t i=0;i<scope.Height();i+=nBinning)
{
// for each row in scope:
// - calculate area for each node in row
// - write row to output file
start=clock();
for (size_t j=0;j<writeScope.Width();j++)
area[j]=WayPoint(nodes+nodesOffset+i*readScope.Width()+j*nBinning).GetAreaDijkstra(limitTime)/areaScaling;
output.UpdateScope(writeScope,area);
writeScope.First(ROW)++;
fprintf(MPI_STDOUT, "\n\trow %u (%u/%u) done in %.2f s\n",i+scope.First(ROW),(i/nBinning)+1,(scope.Height()-1)/nBinning+1,float(clock()-start)/CLOCKS_PER_SEC);
}
output.Close();
delete [] area;
delete [] nodes;
return ERROR_SUCCESS;
}
/*********************************************************************
Copyright 2016, Michael Steinert
This file is part of DrawArea (Dijkstra Routing Algorithm Which
Ascertains REachable Areas).
DrawArea is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
DrawArea is distributed WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DrawArea. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#ifndef A_STAR_INCLUDED
#define A_STAR_INCLUDED
#include <vector>
#include <set>
#include <stdio.h>
#define INVALID_ALTITUDE 0
#define SQRT_2 1.41421356237
// Node is the basic structure of the graph
struct Node {
float altitude; // the altitude at this position
std::vector<struct Edge> neighbors; // vector of edges to all nearest and next-nearest neighbors
Node(size_t alt=0);
void AddNeighbour(int shift, float distance);
};
// Edge is the basic structure to connect nodes of the graph
struct Edge {
static float minTime; // minimum pixel pass time
int shift; // pointer shift to get from node to neighbor
float distance; // temporal distance between the neighboring nodes
Edge(Node* start, int shift, float spacialDistance);
static float Distance(float start, float finish, float spacialDistance); // calculating temporal edge length
};
void CreateGraphEdges(Node* nodes, size_t numCols, size_t numRows);
// MapPoint is basically just a wrapper
// to list, sort and compare WayPoint pointers
struct MapPoint {
class WayPoint* wayPoint;
MapPoint(class WayPoint* wp);
};
bool operator<(const MapPoint& mp1, const MapPoint& mp2);
// WayPoint is a point on the map and stores the
// minimal time to get to the origin
class WayPoint {
private:
Node* node;
float minOriginTime; // currently known shortest time to the origin
public:
WayPoint(Node* node);
~WayPoint();
const WayPoint* GetRoute();
void Expand(class OpenList* openList, class ClosedList* ClosedList, float limitTime);
size_t GetAreaDijkstra(float limitTime);
bool operator>(const WayPoint& wp);
friend bool operator==(const MapPoint& mp, const WayPoint* wp);
};
// all nodes which are partially evaluated
// (a route to the origin is known,
// but not necessarily the shortest)
// are stored in an OpenList
class OpenList {
private:
std::vector<MapPoint> v;
std::vector<MapPoint>::iterator it;
public:
OpenList();
~OpenList();
void Push(WayPoint* wp);
WayPoint* Pop();
void Update();
size_t Size();
bool Empty();
bool Contains(WayPoint* wp);
WayPoint* Select();
};
// all nodes which are fully evaluated
// (the shortest route to the origin is known)
// are stored in a ClosedList
class ClosedList {
private:
std::set<Node*> v;
public:
ClosedList();
~ClosedList();
void Push(Node* node);
size_t Size();
bool Empty();
bool Contains(Node* node);
};
#endif // A_STAR_INCLUDED
/*********************************************************************
Copyright 2016, Michael Steinert
This file is part of DrawArea (Dijkstra Routing Algorithm Which
Ascertains REachable Areas).
DrawArea is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
DrawArea is distributed WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DrawArea. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#ifndef ERROR_H_INCLUDED
#define ERROR_H_INCLUDED
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#define ERROR_SUCCESS 0
#define ERROR_WARNING 0
#ifdef _WIN32
#include <windows.h>
#endif
void InitErrorHandling(const char* appName, const char* heltText, bool consoleErrorOut);
int Error(const char* str, ...);
int Error(int errCode, const char* str, ...);
int Help(const char* msg=NULL);
int Help(int code, const char* msg=NULL);
static const char* HELP_TEXT=NULL;
static const char* APP_NAME=NULL;
static bool CONSOLE_ERROR_OUT=true;
inline void InitErrorHandling(const char* appName, const char* helpText, bool consoleErrorOut)
{
APP_NAME=appName;
HELP_TEXT=helpText;
CONSOLE_ERROR_OUT=consoleErrorOut;
}
// print error message
inline int Error(int code, const char* str, ...)
{
// set prefix to error if error code != 0
va_list pArg;
// get length of error msg to be printed
va_start(pArg,str);
int len=vsnprintf(NULL,0,str,pArg)+1;
va_end(pArg);
// allocate memory and write error message to string
va_start(pArg,str);
char* msg=(char*)malloc(len);
vsprintf(msg,str,pArg);
va_end(pArg);
// output error message
if (CONSOLE_ERROR_OUT)
fprintf(stderr,(code!=ERROR_SUCCESS?"Error: %s":"%s"),msg);
#ifdef _WIN32
else MessageBox(NULL,msg,(code!=0?"Error":"Help"),MB_OK|(code!=0?MB_ICONERROR:MB_ICONINFORMATION));
#endif
free(msg);
return code;
}
// print error message with default error code -1
inline int Error(const char* str, ...)
{
va_list pArg;
va_start(pArg,str);
return Error(-1,str,pArg);
}
// create help message
inline int Help(int code, const char* msg)
{
if (msg)
{
char* msgHelp=(char*)malloc(strlen(HELP_TEXT)+strlen(msg)+2);
strcpy(msgHelp,msg);
strcpy(msgHelp+strlen(msgHelp),"\n\n");
strcpy(msgHelp+strlen(msgHelp),HELP_TEXT);
code=Error(code,msgHelp,APP_NAME);
free(msgHelp);
return code;
}
return Error(code,HELP_TEXT,APP_NAME);
}
inline int Help(const char* msg)
{
return Help(ERROR_SUCCESS,msg);
}
#endif // ERROR_H_INCLUDED
/*********************************************************************
Copyright 2016, Michael Steinert
This file is part of DrawArea (Dijkstra Routing Algorithm Which
Ascertains REachable Areas).
DrawArea is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
DrawArea is distributed WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DrawArea. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#include "DrawArea.h"
// construct Node
Node::Node(size_t alt)
{
altitude=alt;
}
// add a neighboring node and calculate the time to travel to it
// by applying the DistanceFucntion
void Node::AddNeighbour(int shift, float distance)
{
if ((this+shift)->altitude==INVALID_ALTITUDE)
return;
neighbors.push_back(Edge(this,shift,distance));
}
// construct Edge as pointer-shift and temporal distance
float Edge::minTime;
Edge::Edge(Node* start, int shift, float spacialDistance)
{
this->shift=shift;
this->distance=Distance(start->altitude,(start+shift)->altitude,spacialDistance);
}
// calculate travel time along the edge based on altitude difference and spacial distance
float Edge::Distance(float startAltitude, float endAltitude, float spacialDistance)
{
return 0.5*spacialDistance*(startAltitude+endAltitude);
}
// connect each node with its 8 nearest neighbors (if existing)
void CreateGraphEdges(Node* nodes, size_t numCols, size_t numRows)
{
Node* node=nodes;
for (size_t i=0; i<numRows; i++)
for (size_t j=0; j<numCols; j++)
{
node->neighbors.clear();
if (node->altitude==INVALID_ALTITUDE)
{
node++;
continue;
}
if (i>0)
{
node->AddNeighbour(-numCols,1.0); // northern neighbor
if (j>0)
node->AddNeighbour(-numCols-1,SQRT_2); // north-western
if (j+1<numCols)
node->AddNeighbour(-numCols+1,SQRT_2); // north-eastern
}
if (i+1<numRows)
{
node->AddNeighbour(numCols,1.0); // southern neighbor
if (j>0)
node->AddNeighbour(numCols-1,SQRT_2); // south-western
if (j+1<numCols)
node->AddNeighbour(numCols+1,SQRT_2); // south-eastern
}
if (j>0)
node->AddNeighbour(-1,1.0); // western neighbor
if (j+1<numCols)
node->AddNeighbour(1,1.0); // eastern neighbor
node++;
}
}
/*********************************************************************
Copyright 2016, Michael Steinert
This file is part of DrawArea (Dijkstra Routing Algorithm Which
Ascertains REachable Areas).
DrawArea is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
DrawArea is distributed WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DrawArea. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#include "DrawArea.h"
#include <algorithm>
// constructor:
// heap list so we can use heap operations
// to have the MapPoint with the smallest
// estimated distance in front
OpenList::OpenList()
{
make_heap(v.begin(),v.end());
it=v.end();
}
// destructor:
// delete all WayPoints in the list
OpenList::~OpenList()
{
while (!v.empty())
{
delete v.back().wayPoint;
v.pop_back();
}
}
// add MapPoint to list and ensure heap shape
void OpenList::Push(WayPoint* wp)
{
MapPoint mp=new WayPoint(*wp);
v.push_back(mp);
push_heap(v.begin(),v.end());
it=v.end();
}
// remove MapPoint with the smallest
// estimated distance from list and return it
WayPoint* OpenList::Pop()
{
if (v.empty())
return 0;
WayPoint* wp=v.front().wayPoint;
pop_heap(v.begin(),v.end());
v.pop_back();
return wp;
}
// update heap structure of list
void OpenList::Update()
{
make_heap(v.begin(),v.end());
}
// get list size
size_t OpenList::Size()
{
return v.size();
}
// check if list is empty
bool OpenList::Empty()
{
return v.empty();
}
// check if MapPoint is contained in list
bool OpenList::Contains(WayPoint* wp)
{
for (it=v.begin(); it!=v.end(); it++)
if (*it==wp)
return true;
return false;
}
// return a pointer to the WayPoint was wanted
// recently, if it was found
// else return NULL-pointer
WayPoint* OpenList::Select()
{
if (it==v.end())
return NULL;
return it->wayPoint;
}
/*********************************************************************
Copyright 2016, Michael Steinert
This file is part of DrawArea (Dijkstra Routing Algorithm Which
Ascertains REachable Areas).
DrawArea is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
DrawArea is distributed WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DrawArea. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#ifndef PARSECMDLARG_H
#define PARSECMDLARG_H
#include <stdlib.h>
#include <string.h>
/**********************************************************************
This part is general purpose stuff to parse command line arguments.
***********************************************************************/
#ifdef _WIN32
char PATH_SEPARATOR='\\';
#else
char PATH_SEPARATOR='/';
#endif
bool getCmdOption(char** begin, char** end, const char* option);
bool getCmdOption(char** begin, char** end, const char* option, char* value);
bool getCmdOption(char** begin, char** end, const char* option, int& value);
bool getCmdOption(char** begin, char** end, const char* option, float& value);
char* ExtractAppName(const char* cmd);
inline char** find(char** begin, char** end, const char* target)
{
char** itr=begin;
while (itr<end && strcmp(*itr,target)!=0 ) ++itr;
return itr;
}
inline bool getCmdOption(char** begin, char** end, const char* option)
{
return find(begin,end,option)!=end;
}
inline bool getCmdOption(char** begin, char** end, const char* option, char* value)
{
char** itr=find(begin,end,option);
if ( (itr==end) || (++itr==end) ) return false;
strcpy(value,*itr);
return true;
}
inline bool getCmdOption(char** begin, char** end, const char* option, int& value)
{
char str[1024];
char* error;
if (!getCmdOption(begin,end,option,str)) return false;
value=strtol(str,&error,10);
return *error==0;
}
inline bool getCmdOption(char** begin, char** end, const char* option, float& value)
{
char str[1024];
char* error;
if (!getCmdOption(begin,end,option,str)) return false;
value=strtof(str,&error);
return *error==0;
}
inline char* ExtractAppName(char* cmd)
{
char *fileName=cmd;
char *path=strrchr(fileName,PATH_SEPARATOR); // get start pointer of filename w/o path
if (path) fileName=path+1; // set filename start
int cCount=strlen(fileName); // get length of file name
char* ext=strrchr(fileName,'.'); // get start pointer of old extension
if (ext && ext!=fileName) cCount=ext-fileName; // get length of file name w/o extension
char* str=new char[cCount+1]; // create char array to store file name and
strncpy(str,fileName,cCount); // copy old file name w/o extension
str[cCount]='\0';
return str;
}
/**********************************************************************
This part is specially designed for DrawArea.
***********************************************************************/
#include "Error.h"
#include "pnm.h"
struct CmdArgStruct
{
Scope scope;
float areaScaling;
int nBinning;
int nColumns;
float limitTime;
float minTime;
bool noCalc;
int nRows;
int nThreads;
char inputName[256];
char outputName[256];
int Evaluate(int argc, char** argv);
};
// evaluated command line arguments
inline int CmdArgStruct::Evaluate(int argc, char** argv)
{
// set input file
if (argc<2)
return Help(-1,"Missing input file.");
strcpy(inputName,argv[argc-1]);
// set area scaling
areaScaling=1.0;
if (getCmdOption(argv+1,argv+argc,"-a",areaScaling) || getCmdOption(argv+1,argv+argc,"--area",areaScaling))
if (!areaScaling>0)
return Help(-1,"Area scaling factor must be a positive.");
// set binning
nBinning=1;
if (getCmdOption(argv+1,argv+argc,"-b",nBinning) || getCmdOption(argv+1,argv+argc,"--binning",nBinning))
if (nBinning<1)
return Help(-1,"Binning must be a positive.");
areaScaling*=nBinning*nBinning;
// set number of columns to split into
nColumns=0;
if (getCmdOption(argv+1,argv+argc,"-c") || getCmdOption(argv+1,argv+argc,"--cols"))
{
nColumns=-1;
if (getCmdOption(argv+1,argv+argc,"-c",nColumns) || getCmdOption(argv+1,argv+argc,"--cols",nColumns))
if (nColumns<1)
return Help(-1,"Number of columns must be a positive.");
}
// set time limit
if (getCmdOption(argv+1,argv+argc,"-l",limitTime) || getCmdOption(argv+1,argv+argc,"--limit",limitTime))
{
if (!limitTime>0)
return Help(-1,"Travel time limit must be a positive.");
}
else return Help(-1,"Missing travel time limit.");
// set minimum pixel pass time
if (getCmdOption(argv+1,argv+argc,"-m",minTime) || getCmdOption(argv+1,argv+argc,"--min-time",minTime))
{
if (!minTime>0)
return Help(-1,"Minimum time to pass single pixel must be a positive.");
}
else return Help(-1,"Missing minimum pixel time.");
// user requested only scope splitting but no map calculation
noCalc=false;
if (getCmdOption(argv+1,argv+argc,"-n") || getCmdOption(argv+1,argv+argc,"--no-calc"))
noCalc=true;
// set output name
if (!getCmdOption(argv+1,argv+argc,"-o",outputName) && !getCmdOption(argv+1,argv+argc,"--output",outputName))
{
strcpy(outputName,"area_of_\0");
strcpy(outputName+strlen(outputName),inputName);
}
// set number of rows to split into
nRows=0;
if (getCmdOption(argv+1,argv+argc,"-r") || getCmdOption(argv+1,argv+argc,"--rows"))
{
nRows=-1;
if (getCmdOption(argv+1,argv+argc,"-r",nRows) || getCmdOption(argv+1,argv+argc,"--cols",nRows))
if (nRows<1)
return Help(-1,"Number of rows must be a positive.");
}
// set user defined scope
scope=Scope(-1,-1);
char strScope[256];
if (getCmdOption(argv+1,argv+argc,"-s",strScope) || getCmdOption(argv+1,argv+argc,"--scope",strScope))
if (sscanf(strScope,"%ldx%ld+%ld+%ld",&scope.Width(),&scope.Height(),&scope.First(COLUMN),&scope.First(ROW))!=4)
return Help(-1,"Scope must be of format INTxINT+INT+INT.");
// set number of threads
nThreads=1;
if (getCmdOption(argv+1,argv+argc,"-t",nThreads) || getCmdOption(argv+1,argv+argc,"--threads",nThreads))
if (nThreads<1)
return Help(-1,"Number of threads must be a positive.");
// calculate number of rows and columns to split into
if (nRows>0 && nColumns>0)
nThreads=nRows*nColumns;
else
{
if (nColumns)
{
nRows=1;
nThreads=nColumns=(nColumns>0?nColumns:nThreads);
}
else
{
nColumns=1;
nThreads=nRows=(nRows>0?nRows:nThreads);
}
}
return ERROR_SUCCESS;
}
#endif
/*********************************************************************
Copyright 2016, Michael Steinert
This file is part of DrawArea (Dijkstra Routing Algorithm Which
Ascertains REachable Areas).
DrawArea is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
DrawArea is distributed WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DrawArea. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#include "DrawArea.h"
#include <algorithm>
// constructor:
// connect the MapPoint with a WayPoint
MapPoint::MapPoint(WayPoint* wp)
{
wayPoint=wp;
}
// compare two MapPoints
// this is needed to heap the open list
bool operator<(const MapPoint& mp1, const MapPoint& mp2)
{
return *(mp1.wayPoint)>*(mp2.wayPoint);
}
// constructor:
// connect the WayPoint with a node
//
// this is quite useful because with this
// you can use a Node pointer in functions
// that require a WayPoint
WayPoint::WayPoint(Node* node)
{
this->node=node;
minOriginTime=0.0;
}
// destructor (see constructor ClosedList)
WayPoint::~WayPoint() {}
// evaluate all neighboring nodes by estimating time to travel through this nodes
void WayPoint::Expand(OpenList* openList, ClosedList* closedList, float limitTime)
{
// loop through all neighboring nodes
for (std::vector<Edge>::iterator it=node->neighbors.begin(); it!=node->neighbors.end(); ++it)
{
// if current neighbor is fully evaluated (contained in closedList)
// continue with the next neighbor
if (closedList->Contains(node+it->shift))
continue;
WayPoint currNeighbour(node+it->shift);
WayPoint* pCurrNeighbour=&currNeighbour;
// check if node is already contained in openList
// if so, use WayPoint in openList
bool inOpenList=openList->Contains(pCurrNeighbour);
if (inOpenList)
pCurrNeighbour=openList->Select();
// calculate time to origin as sum of previous node time and
// time to get from there to this node
float originTime=minOriginTime+it->distance;
// if time limit is exceeded continue with next neighbor
if (originTime>limitTime)
continue;
// if shorter route is known continue with next neighbor
if (inOpenList && originTime>=pCurrNeighbour->minOriginTime)
continue;
// refresh times and shortest route previous node
pCurrNeighbour->minOriginTime=originTime;
// update openList
if (inOpenList)
openList->Update();
else
openList->Push(pCurrNeighbour);
}
}
// this function is looking for all nodes which could be reached in time limit
// by using the Dijkstra algorithm
size_t WayPoint::GetAreaDijkstra(float limitTime)
{
if (node->altitude==INVALID_ALTITUDE)
return 0;
minOriginTime=0.0;
// define empty closed and open list
ClosedList closedList;
OpenList openList;
// start at origin
openList.Push(this);
// loop until all nodes that could be reached in time limit are fully evaluated
while (!openList.Empty())
{
// use the WayPoint with the lowest travel time
WayPoint* top=openList.Pop();
// add the WayPoint to the closed list to mark it as fully evaluated
closedList.Push(top->node);
// evaluate all neighboring nodes
top->Expand(&openList,&closedList,limitTime);
delete top;
}
return closedList.Size();
}
bool WayPoint::operator>(const WayPoint& wp)
{
return minOriginTime>wp.minOriginTime;
}
// check if MapPoint node and a Node pointer pointing to the same Node
bool operator==(const MapPoint& mp, const WayPoint* wp)
{
return wp->node==mp.wayPoint->node;
}
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment