ngmf.util.cosu.Fast Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of oms Show documentation
Show all versions of oms Show documentation
Object Modeling System (OMS) is a pure Java object-oriented framework.
OMS v3.+ is a highly interoperable and lightweight modeling framework for component-based model and simulation development on multiple platforms.
/*
* $Id: Fast.java 50798ee5e25c 2013-01-09 [email protected] $
*
* This file is part of the Object Modeling System (OMS),
* 2007-2012, Olaf David and others, Colorado State University.
*
* OMS is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, version 2.1.
*
* OMS is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with OMS. If not, see .
*/
package ngmf.util.cosu;
import java.util.Arrays;
/**
*
* @author joy,od
*/
public class Fast {
public static void main(String[] args) {
// FAST operation and external model parameters
int N = 2500; // number of samples
double M = 4; // number of terms in the partial variances summation
int npar = 4; // number of parameters
double LB[] = {1, 0, 0, 0};
double UB[] = {2, 1, 1, 1};
// Add warning if N is less than minimum required realizations
if (N < 2000) {
throw new RuntimeException("The number of realizations is less than 2000.");
}
double wi = Math.floor(N / (2 * M));
double m2 = Math.floor(wi / (2 * M));
double r = Math.floor((m2) / npar);
double[] w = new double[npar];
if (r < 1) {
for (int i = 0; i <= npar - 1; i++) {
w[i] = 1;
}
} else {
double t = Math.floor(m2 / npar);
w[0] = 1;
for (int i = 1; i <= npar - 1; i++) {
w[i] = 1 + i * t;
}
}
int k1 = 0;
double[][] w2 = new double[npar][w.length];
for (int i = 0; i <= npar - 1; i++) {
for (int j = 0; j <= w.length - 1; j++) {
if (j == k1) {
w2[i][j] = wi;
} else {
w2[i][j] = w[j];
}
}
k1++;
}
double inc = 2 * Math.PI / N;
double[] s = new double[N];
s[0] = -Math.PI;
for (int i = 1; i < s.length; i++) {
s[i] = s[i - 1] + inc;
}
double[][] x = new double[N][npar];
double[] y = new double[N];
double[] V = new double[npar];
double[] VT = new double[npar];
double[] Ak = new double[(int) Math.floor((N - 1) / 2)];
double[] Bk = new double[(int) Math.floor((N - 1) / 2)];
double[] S_par = new double[npar];
double[] Vex = new double[npar];
double[] Sex_par = new double[npar];
// System.out.println((int) Math.floor((N - 1) / 2));
// System.out.println((N - 1) / 2);
for (int h = 0; h <= npar - 1; h++) {
// Compute realizations
for (int j = 0; j <= N - 1; j++) {
for (int i = 0; i <= npar - 1; i++) {
double p = 0.5 + Math.asin(Math.sin(w2[h][i] * s[j])) / Math.PI;
p = (double) Math.round(p * 10000) / 10000;
x[j][i] = p * (UB[i] - LB[i]) + LB[i];
}
y[j] = run_model(x[j][0], x[j][1], x[j][2], x[j][3]);
}
// Compute total variance
V[h] = 0;
for (int k = 1; k <= ((N - 1) / 2); k++) {
double A = 0, B = 0;
for (int j = 0; j <= N - 1; j++) {
if (j == 0) {
A = y[j] * Math.cos(s[j] * k);
B = y[j] * Math.sin(s[j] * k);
} else {
A += y[j] * Math.cos(s[j] * k); //#ok
B += y[j] * Math.sin(s[j] * k); //#ok
}
A = (double) Math.round(A * 10000) / 10000;
B = (double) Math.round(B * 10000) / 10000;
}
Ak[k - 1] = A * 2 / N;
Bk[k - 1] = B * 2 / N;
V[h] = V[h] + Math.pow((A * 2 / N), 2) + Math.pow((B * 2 / N), 2);
V[h] = (double) Math.round(V[h] * 10000) / 10000;
}
VT[h] = V[h] / 2;
//Compute partial variance
V[h] = 0;
for (int q = 1; q <= M; q++) {
V[h] = V[h] + Math.pow(Ak[(int) (q * w2[h][h]) - 1], 2)
+ Math.pow(Bk[(int) (q * w2[h][h]) - 1], 2);
V[h] = (double) Math.round(V[h] * 10000) / 10000;
}
V[h] = V[h] / 2;
S_par[h] = V[h] / VT[h];
S_par[h] = (double) Math.round(S_par[h] * 1000000) / 1000000;
//Compute Extended partial variance
Vex[h] = 0;
for (int q = 1; q <= M; q++) {
for (int c = 0; c <= npar - 1; c++) {
if (c != h) {
Vex[h] = Vex[h] + Math.pow(Ak[(int) (q * w2[h][c]) - 1], 2)
+ Math.pow(Bk[(int) (q * w2[h][c]) - 1], 2);
Vex[h] = (double) Math.round(Vex[h] * 10000) / 10000;
}
}
}
Vex[h] /= 2;
Sex_par[h] = (1 - Vex[h] / VT[h]);
Sex_par[h] = (double) Math.round(Sex_par[h] * 1000000) / 1000000;
}
//print out S values
System.out.println("S =" + Arrays.toString(S_par));
System.out.println("ST =" + Arrays.toString(Sex_par));
}
static double run_model(double x1, double x2, double x3, double x4) {
double y = x1 * x2 + x3 + x4;
return y;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy