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 jgt-oms3 Show documentation
There is a newer version: 0.8.1
/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
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;
    }
}