use of org.hortonmachine.gears.utils.sorting.QuickSortAlgorithm in project hortonmachine by TheHortonMachine.
the class OmsMultiTca method process.
@Execute
public void process() {
if (!concatOr(outMultiTca == null, doReset)) {
return;
}
checkNull(inPit, inFlow, inCp9);
HashMap<String, Double> regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inFlow);
int cols = regionMap.get(CoverageUtilities.COLS).intValue();
int rows = regionMap.get(CoverageUtilities.ROWS).intValue();
// pm.message();
@SuppressWarnings("unused") int ipos, jpos, i, j, ncicli = 0;
double sum, delta, pos;
// create new matrix
double[] elevationArray = new double[cols * rows];
double[] indexOfElevation = new double[cols * rows];
// pm.message();
RandomIter flowIter = CoverageUtilities.getRandomIterator(inFlow);
RandomIter pitIter = CoverageUtilities.getRandomIterator(inPit);
RandomIter cp9Iter = CoverageUtilities.getRandomIterator(inCp9);
WritableRaster alreadyDonePixelWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, 0.0);
WritableRandomIter alreadyDoneIter = RandomIterFactory.createWritable(alreadyDonePixelWR, null);
WritableRaster multiTcaWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, 1.0);
WritableRandomIter multiTcaIter = RandomIterFactory.createWritable(multiTcaWR, null);
/*
* store the value of elevation in an array
*/
for (int t = 0; t < rows; t++) {
for (int s = 0; s < cols; s++) {
elevationArray[((t) * cols) + s] = pitIter.getSampleDouble(s, t, 0);
indexOfElevation[((t) * cols) + s] = ((t) * cols) + s + 1;
}
}
// pm.message();
try {
QuickSortAlgorithm sortAlgorithm = new QuickSortAlgorithm(pm);
sortAlgorithm.sort(elevationArray, indexOfElevation);
} catch (Exception e) {
e.printStackTrace();
}
/*
* start to working with the highest value of elevation.
*/
for (int l = cols * rows - 1; l >= 0; l--) {
ncicli = cols * rows - l;
if (elevationArray[l] <= 0) {
break;
} else {
pos = indexOfElevation[l];
// extract the index of the matrix from the arrays index.
i = (int) (pos - 1) % cols;
j = (int) (pos - 1) / cols;
if (alreadyDoneIter.getSampleDouble(i, j, 0) == 0.0) {
alreadyDoneIter.setSample(i, j, 0, 1.0);
if (cp9Iter.getSampleDouble(i, j, 0) == 10 || cp9Iter.getSampleDouble(i, j, 0) == 20 || cp9Iter.getSampleDouble(i, j, 0) == 30 || cp9Iter.getSampleDouble(i, j, 0) == 40 || cp9Iter.getSampleDouble(i, j, 0) == 50 || cp9Iter.getSampleDouble(i, j, 0) == 60) {
sum = 0;
for (int k = 1; k <= 8; k++) {
ipos = i + dir[k][0];
jpos = j + dir[k][1];
delta = pitIter.getSampleDouble(i, j, 0) - pitIter.getSampleDouble(ipos, jpos, 0);
if (delta == 0) {
if (alreadyDoneIter.getSampleDouble(ipos, jpos, 0) == 0.0 && flowIter.getSampleDouble(ipos, jpos, 0) == dirIn[k][2]) {
resolveFlat(ipos, jpos, cols, rows, pitIter, multiTcaIter, alreadyDoneIter, flowIter, cp9Iter);
}
}
if (delta > 0.0 && pitIter.getSampleDouble(ipos, jpos, 0) > 0.0) {
sum += delta;
}
}
for (int k = 1; k <= 8; k++) {
ipos = i + dir[k][0];
jpos = j + dir[k][1];
delta = pitIter.getSampleDouble(i, j, 0) - pitIter.getSampleDouble(ipos, jpos, 0);
if (delta > 0.0 && pitIter.getSampleDouble(ipos, jpos, 0) > 0.0) {
multiTcaIter.setSample(ipos, jpos, 0, multiTcaIter.getSampleDouble(ipos, jpos, 0) + multiTcaIter.getSampleDouble(i, j, 0) * (delta / sum));
} else if (delta == 0.0 && flowIter.getSampleDouble(i, j, 0) == dirIn[k][2]) {
multiTcaIter.setSample(ipos, jpos, 0, multiTcaIter.getSampleDouble(ipos, jpos, 0) + multiTcaIter.getSampleDouble(i, j, 0));
}
}
} else if (cp9Iter.getSampleDouble(i, j, 0) == 70 || cp9Iter.getSampleDouble(i, j, 0) == 80 || cp9Iter.getSampleDouble(i, j, 0) == 90) {
for (int k = 1; k <= 8; k++) {
ipos = i + dir[k][0];
jpos = j + dir[k][1];
double delta2 = pitIter.getSampleDouble(i, j, 0) - pitIter.getSampleDouble(ipos, jpos, 0);
if (delta2 == 0) {
if (alreadyDoneIter.getSampleDouble(ipos, jpos, 0) == 0.0 && flowIter.getSampleDouble(ipos, jpos, 0) == dirIn[k][2]) {
resolveFlat(ipos, jpos, cols, rows, pitIter, multiTcaIter, alreadyDoneIter, flowIter, cp9Iter);
}
}
}
for (int k = 1; k <= 8; k++) {
ipos = i + dir[k][0];
jpos = j + dir[k][1];
if (flowIter.getSampleDouble(i, j, 0) != 10 && flowIter.getSampleDouble(i, j, 0) == dir[k][2]) {
multiTcaIter.setSample(ipos, jpos, 0, multiTcaIter.getSampleDouble(ipos, jpos, 0) + multiTcaIter.getSampleDouble(i, j, 0));
break;
}
}
}
}
}
}
for (int t = 0; t < rows; t++) {
for (int s = 0; s < cols; s++) {
if (isNovalue(cp9Iter.getSampleDouble(s, t, 0)) || isNovalue(flowIter.getSampleDouble(s, t, 0)))
multiTcaIter.setSample(s, t, 0, HMConstants.doubleNovalue);
}
}
outMultiTca = CoverageUtilities.buildCoverage("multiTca", multiTcaWR, regionMap, inFlow.getCoordinateReferenceSystem());
}
use of org.hortonmachine.gears.utils.sorting.QuickSortAlgorithm in project hortonmachine by TheHortonMachine.
the class GeometryUtilities method sortHorizontal.
public static void sortHorizontal(Coordinate[] coordinates) {
QuickSortAlgorithm sorter = new QuickSortAlgorithm(new DummyProgressMonitor());
double[] x = new double[coordinates.length];
double[] y = new double[coordinates.length];
for (int i = 0; i < x.length; i++) {
x[i] = coordinates[i].x;
y[i] = coordinates[i].y;
}
sorter.sort(x, y);
for (int i = 0; i < x.length; i++) {
coordinates[i].x = x[i];
coordinates[i].y = y[i];
}
}
use of org.hortonmachine.gears.utils.sorting.QuickSortAlgorithm in project hortonmachine by TheHortonMachine.
the class CoupledFieldsMoments method process.
public double[][] process(RenderedImage map1RI, RenderedImage map2RI, int pBins, int pFirst, int pLast, IHMProgressMonitor pm, int binmode) {
if (map2RI == null) {
map2RI = map1RI;
}
GearsMessageHandler msg = GearsMessageHandler.getInstance();
pm.message(msg.message("cb.vectorize"));
double[] U = vectorizeDoubleMatrix(map1RI);
double[] T = null;
QuickSortAlgorithm t = new QuickSortAlgorithm(pm);
if (map2RI == null) {
T = U;
t.sort(U, (double[]) null);
} else {
T = vectorizeDoubleMatrix(map2RI);
t.sort(U, T);
}
SplitVectors theSplit = new SplitVectors();
int num_max = 1000;
/*
* if (bintype == 1) {
*/
pm.message(msg.message("cb.splitvector"));
split2realvectors(U, T, theSplit, pBins, num_max, pm);
/*
* } else { delta = FluidUtils.exponentialsplit2realvectors(U, T,
* theSplit, N, num_max, base); }
*/
pm.message(msg.message("cb.creatematrix"));
double[][] outCb = new double[theSplit.splitIndex.length][pLast - pFirst + 3];
// kept for future expansion
binmode = 1;
if (// always true for now, other modes not implemented yet
binmode == 1) {
for (int h = 0; h < theSplit.splitIndex.length; h++) {
outCb[h][0] = calculateNthMoment(theSplit.splitValues1[h], (int) theSplit.splitIndex[h], 0.0, 1.0, pm);
outCb[h][1] = theSplit.splitIndex[h];
outCb[h][2] = calculateNthMoment(theSplit.splitValues2[h], (int) theSplit.splitIndex[h], 0.0, 1.0, pm);
if (pFirst == 1)
pFirst++;
for (int k = pFirst; k <= pLast; k++) {
outCb[h][k - pFirst + 3] = calculateNthMoment(theSplit.splitValues2[h], (int) theSplit.splitIndex[h], outCb[h][1], (double) k, pm);
}
}
}
return outCb;
}
use of org.hortonmachine.gears.utils.sorting.QuickSortAlgorithm in project hortonmachine by TheHortonMachine.
the class OmsJami method extractFromStationFeatures.
/**
* Fills the elevation and id arrays for the stations, ordering in ascending
* elevation order.
*
* @throws Exception
* in the case the sorting gives problems.
*/
private void extractFromStationFeatures() throws Exception {
int stationIdIndex = -1;
int stationElevIndex = -1;
pm.beginTask("Filling the elevation and id arrays for the stations, ordering them in ascending elevation order.", stationCoordinates.size());
for (int i = 0; i < stationCoordinates.size(); i++) {
pm.worked(1);
SimpleFeature stationF = stationFeatures.get(i);
Coordinate stationCoord = stationCoordinates.get(i);
if (stationIdIndex == -1) {
SimpleFeatureType featureType = stationF.getFeatureType();
stationIdIndex = featureType.indexOf(fStationid);
stationElevIndex = featureType.indexOf(fStationelev);
if (stationIdIndex == -1) {
throw new IllegalArgumentException("Could not find the field: " + fStationid);
}
if (stationElevIndex == -1) {
throw new IllegalArgumentException("Could not find the field: " + fStationelev);
}
}
int id = ((Number) stationF.getAttribute(stationIdIndex)).intValue();
double elev = ((Number) stationF.getAttribute(stationElevIndex)).doubleValue();
statElev[i] = elev;
statId[i] = id;
stationId2CoordinateMap.put(id, stationCoord);
}
pm.done();
// sort
QuickSortAlgorithm qsA = new QuickSortAlgorithm(pm);
qsA.sort(statElev, statId);
}
use of org.hortonmachine.gears.utils.sorting.QuickSortAlgorithm in project hortonmachine by TheHortonMachine.
the class NetworkBuilder method geoSewer.
/**
* Project the net.
*
* <p>
* it calculate some properties of each pipe, for instance the diameter, the
* depth, the discharge.
* </p>
* <p>
* The Network is a collection of {@link Pipe}, and the discharge is
* evalutate using a pliviometric curve.
* </p>
* <p>
* There is two loops:
* <ol>
* <li>The first on the head pipes.
* <li>The second on the pipes which are internal.
* </ol>
* </p>
*/
@Override
public void geoSewer() throws Exception {
/* l Tratto che si sta progettando. */
int l;
/*
* Estensione del sottobacino con chiusura nel tratto l che si sta
* progettando.
*/
double totalarea;
/*
* Passo con cui variare t e tp, nella ricerca della portata massima di
* progetto.
*/
double dtp;
/*
* Tiene traccia dei diametri utilizzati e fa in modo che procedendo
* verso valle non vi siano restringimenti.
*/
double maxd;
/* r di tentativo, dove r = tp / k. */
double oldr;
/* r Valore finale di r ( tp / k ). */
double r = 0;
/* No Vale L / ( k * c ). */
double No = 0;
/*
* Estremo inferiore dell'intervallo che contiene la radice ricercata,
* ossia la r* che massimizza la portata.
*/
double inf;
/*
* Estremo superiore dell'intervallo che contiene la r*, che massimizza
* la portata.
*/
double sup;
/* Portata di tentativo. */
double oldQ;
/*
* matrice che per ciascun area non di testa, contiene i dati geometrici
* degli stati a monte, che direttamente o indirettamente, drenano in
* esso
*/
double[][] net;
/*
* contiene la magnitude dei vari stati.
*/
double[] magnitude = new double[networkPipes.length];
/*
* vettore che contiene l'indice dei versanti.
*/
double[] one = new double[networkPipes.length];
/*
* vettore che contiene gli stati riceventi, compresa almeno un'uscita
*/
double[] two = new double[networkPipes.length];
for (int i = 0; i < networkPipes.length; i++) {
/* Indice degli stati */
one[i] = i;
/* Indice degli stati riceventi, compresa almeno un'uscita */
two[i] = networkPipes[i].getIndexPipeWhereDrain();
}
/* Calcola la magnitude di ciascun stato */
Utility.pipeMagnitude(magnitude, two, pm);
/*
* Calcola la magnitude di
* ciascun stato
*/
double tolerance = networkPipes[0].getTolerance();
/* al vettore two vengono assegnati gli elementi di magnitude */
for (int i = 0; i < two.length; i++) {
two[i] = magnitude[i];
}
/*
* Ordina gli elementi del vettore magnitude in ordine crescente, e
* posiziona nello stesso ordine gli elementi di one
*/
QuickSortAlgorithm t = new QuickSortAlgorithm(pm);
t.sort(magnitude, one);
/* ----- INIZIO DIMENSIONAMENTO DELLE AREE DI TESTA ----- */
/*
* qup Indice del tratto a cui corrisponde il diametro massimo, quando
* si analizza un sottobacino.
*/
int[] qup = new int[1];
int k = 0;
// vettore che contiene i ritardi locali dell'onda.
double[] localdelay = new double[magnitude.length];
// tratto che si sta analizzando o progettando
l = (int) one[k];
pm.beginTask(msg.message("trentoP.begin"), networkPipes.length - 1);
int jMax = networkPipes[0].getjMax();
double c = networkPipes[0].getC();
double g = networkPipes[0].getG();
double tau = networkPipes[0].getTau();
minDischarge = networkPipes[0].getMinDischarge();
while (magnitude[k] == 1) {
/*
* Serve per tener conto della forma, piu o meno allungata, delle
* aree drenanti
*/
if (networkPipes[l].getAverageResidenceTime() >= 0) {
/*
* La formula 1.7 ( k = alfa * S ^ beta / ( ksi ^ b * s ^ GAMMA
* ) k tempo di residenza [ min ]
*/
// TODO: check formulas!!
networkPipes[l].residenceTime = ((HOUR2MIN * networkPipes[l].getAverageResidenceTime() * pow(networkPipes[l].getDrainArea() / METER2CM, exponent)) / (pow(networkPipes[l].getRunoffCoefficient(), esp) * pow(networkPipes[l].getAverageSlope(), gamma)));
} else {
/*
* Considero solo l 'acqua che drena dalla strada k = alfa * S /
* L * i ^ GAMMA [ min ]
*/
networkPipes[l].residenceTime = (-networkPipes[l].getDrainArea() * HOUR2MIN * networkPipes[l].getAverageResidenceTime() / networkPipes[l].getLength());
}
maxd = 0;
/*
* Velocita media di primo tentativo nel tratto da progettare [m/s]
*/
networkPipes[l].meanSpeed = (1.0);
// r di primo tentativo [adimensional]
oldr = 1.0;
/*
* Estremo inferiore da adottare nella ricerca della r* che
* massimizza la portata.
*/
inf = 0.1;
/*
* ----- INIZIO CICLO FOR PER LA PROGETTAZIONE DEL TRATTO PARTENDO
* DA UNA r e celerita DI PRIMO TENTATIVO -----
*/
int j;
for (j = 0; j <= jMax; j++) {
/*
* L / ku Calcolato in funzione della velocita di primo
* tentativo . No sara ricalcolato finche non si avra una
* convergenza di r .
*/
No = networkPipes[l].getLength() / (MINUTE2SEC * networkPipes[l].residenceTime * celerityfactor * networkPipes[l].meanSpeed);
// Estremo superiore da adottare nella ricerca della r*.
sup = 2 * No + 5;
/*
* tp* [min] che da origine alla massima portata, calcolato come
* r*k
*/
R_F rfFunction = new R_F();
rfFunction.setParameters(n, No);
r = RootFindingFunctions.bisectionRootFinding(rfFunction, inf, sup, tolerance, jMax, pm);
networkPipes[l].tP = (r * networkPipes[l].residenceTime);
/*
* coefficiente udometrico calcolato con la formula 2.17 u [ l /
* s * ha ] o l/min/ha???
*/
networkPipes[l].coeffUdometrico = (networkPipes[l].getRunoffCoefficient() * a * pow(networkPipes[l].tP, n - 1) * (1 + MINUTE2SEC * celerityfactor * networkPipes[l].meanSpeed * networkPipes[l].tP / networkPipes[l].getLength() - 1 / No * log(exp(No) + exp(r) - 1)) * 166.6666667);
/*
* Portata Q [ l / s ]
*/
// TODO: changed here!!
networkPipes[l].discharge = (networkPipes[l].coeffUdometrico * networkPipes[l].getDrainArea());
networkPipes[l].designPipe(diameters, tau, g, maxd, c, strBuilder);
/*
* La r e stata determinata per via iterativa con la precisione
* richiesta , allora esce dal ciclo for
*/
if (abs((r - oldr) / oldr) <= tolerance) {
/*
* la j che tiene conto del numero di iterazioni viene
* settata a 0
*/
j = 0;
break;
}
/*
* non si e arrivati alla convergenza di r, quindi si usa la
* nuova r per un' ulteriore iterazione
*/
oldr = r;
}
/*
* Si e usciti dal precedente ciclo for perche si e superato il
* numero massimo di iterazioni JMAX ammesse senza arrivare alla
* convergenza di
*/
if (j != 0) {
pm.errorMessage(msg.message("trentoP.error.conv"));
}
/*
* t * [ min ] tempo in cui si verifica la massima tra le portata
* massime all 'uscita del tratta appena progettato
*/
networkPipes[l].tQmax = (networkPipes[l].residenceTime * log(exp(No) + exp(r) - 1));
/*
* L / u [ min ] ritardo locale dell 'onda di piena
*/
localdelay[l] = (networkPipes[l].getLength()) / (celerityfactor * MINUTE2SEC * networkPipes[l].meanSpeed);
// Ac [ha] superfice servita
networkPipes[l].totalSubNetArea = networkPipes[l].getDrainArea();
// Mean length of upstream net [m] (=length of pipe)
networkPipes[l].totalSubNetLength = networkPipes[l].getLength();
networkPipes[l].meanLengthSubNet = networkPipes[l].getLength();
// Passo allo stato successivo
k++;
if (k < magnitude.length) {
/*
* Il prossimo tratto da progettare, ovviamente se avra
* magnitude=1
*/
l = (int) one[k];
} else {
break;
}
pm.worked(1);
}
/*
* ----- FINE CICLO WHILE PER LA PROGETTAZIONE DELLE AREE DI TESTA
*/
dtp = tDTp;
if (dtp > 0.5) {
dtp = 0.5;
/*
* passo temporale con cui valutare le portate quando si
* ricerca la portata massima di progetto [min]
*/
strBuilder.append(msg.message("trentoP.warning.timestep"));
}
/*
* ----- INIZIO CICLO WHILE PER LA PROGETTAZIONE DELLE AREE NON DI TESTA
* -----
*
* Magnitude > 1 AREE NON DI TESTA
*/
while (k < magnitude.length) {
/*
* Crea una matrice net[k][9], dove k e pari al numero di stati,
* che direttamente o indirettamente , drenano nello stato
*/
net = new double[(int) (magnitude[k] - 1)][9];
/*
* Serve per tener conto della forma, piu o meno allungata, delle
* aree drenanti
*/
if (networkPipes[l].getAverageResidenceTime() >= 0) {
/*
* La formula 1.7 ( k = alfa * S ^ beta * i ^ GAMMA ) k tempo di
* residenza [ min ]
*/
networkPipes[l].residenceTime = ((HOUR2MIN * networkPipes[l].getAverageResidenceTime() * pow(networkPipes[l].getDrainArea() / METER2CM, exponent)) / (pow(networkPipes[l].getRunoffCoefficient(), esp) * pow(networkPipes[l].getAverageSlope(), gamma)));
} else {
// k tempo di residenza [ min ]
networkPipes[l].residenceTime = (-networkPipes[l].getDrainArea() * networkPipes[l].getAverageResidenceTime() / networkPipes[l].getLength());
}
// Restituisce l'area del subnetwork che si chiude in l
totalarea = scanNetwork(k, l, one, net, qup);
// Diametro massimo riscontrato nel subnetwork analizzato
maxd = networkPipes[qup[0]].diameter;
/*
* Velocita media di primo tentativo [m/s]
*/
networkPipes[l].meanSpeed = (1.0);
/*
* Portata di primo tentativp [l/s]
*/
networkPipes[l].discharge = (minDischarge);
/*
* ----- INIZIO CICLO DO WHILE (progettare fino alla convergenza
* della Q) -----
*/
do {
oldQ = networkPipes[l].discharge;
networkPipes[l].discharge = (0);
/*
* L / u [ min ]
*/
localdelay[l] = networkPipes[l].getLength() / (celerityfactor * MINUTE2SEC * networkPipes[l].meanSpeed);
/*
* Aggiorna i ritardi nella matrice net, includendo il ritardo
* relativo allo stato che si sta progettando. Questo perche la
* celerita nell'ultimo tratto non e nota a priori, ma verra
* calcolata iteraivamente.
*/
for (int i = 0; i < net.length; i++) {
net[i][2] += localdelay[l];
}
/*
* Restituisce il contributo alla portata dello stato che si sta
* progettando
*/
discharge(l, dtp, net, localdelay);
/*
* Risetta i ritardi, toglieno il ritardo relativo allo stato
* che si sta progettando
*/
for (int i = 0; i < net.length; i++) {
net[i][2] -= localdelay[l];
}
networkPipes[l].designPipe(diameters, tau, g, maxd, c, strBuilder);
/*
* finche' si arriva alla convergenza della portata Q
*/
} while (abs(oldQ - networkPipes[l].discharge) / oldQ > epsilon);
/*
* Coefficiente udometrico u [l/(s*ha )]
*/
networkPipes[l].coeffUdometrico = (networkPipes[l].discharge / totalarea);
/* Ac [ha] */
networkPipes[l].totalSubNetArea = totalarea;
// Mean length of upstream net [ m ]
networkPipes[l].meanLengthSubNet = ModelsEngine.meanDoublematrixColumn(net, 1);
/*
* Variance of lengths of upstream net [ m ^ 2 ]
*/
networkPipes[l].varianceLengthSubNet = ModelsEngine.varianceDoublematrixColumn(net, 1, networkPipes[l].meanLengthSubNet);
/* Passo allo stato successivo */
k++;
/* se non sono arrivato alla fine */
if (k < magnitude.length) {
/* Prossimo stato da progettare */
l = (int) one[k];
} else {
break;
}
pm.worked(1);
}
resetDepths(two);
}
Aggregations