Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Untitled

By: a guest on May 14th, 2012  |  syntax: None  |  size: 9.01 KB  |  hits: 15  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. package KKProgramming;
  2.  
  3. import java.util.ArrayList;
  4.  
  5. import Jama.Matrix;
  6.  
  7. public class Obliczenia {
  8.         private double r;
  9.         private double mi;
  10.         private double epsilon;
  11.         private double alfa;
  12.         private double t0;
  13.         private double gamma;
  14.         private double beta;
  15.  
  16.         private ArrayList<Zasob> supplies;
  17.         private ArrayList<Matrix> buildings;
  18.         public final static int BACKTRACKING = 1;
  19.         public final static int GOLDENRATIO = 2;
  20.         public final static int HOOK = 3;
  21.  
  22.         // statystyki
  23.         public long hookAll = 0;
  24.         public long hookWork = 0;
  25.         public double hookResult = 0.;
  26.  
  27.         public int interiorAll = 0;
  28.         public double newtonAv = 0.;
  29.         public int newtonAll = 0;
  30.         public double interiorResult = 0;
  31.  
  32.         public void clearBuildings() {
  33.                 buildings.clear();
  34.         }
  35.  
  36.         public ArrayList<Matrix> getBuild() {
  37.                 return buildings;
  38.         }
  39.  
  40.         public void setVariables(double r, double t, double gamma, double epsilon,
  41.                         double mi, double alpha) {
  42.                 this.r = r;
  43.                 this.t0 = t;
  44.                 this.gamma = gamma;
  45.                 this.epsilon = epsilon;
  46.                 this.mi = mi;
  47.                 this.alfa = alpha;
  48.         }
  49.  
  50.         private Matrix arrayToMatrix(ArrayList<Matrix> array) {
  51.                 int n = array.size();
  52.                 Matrix m = new Matrix(n * 2, 1);
  53.                 int i = 0;
  54.                 for (Matrix a : array) {
  55.                         m.set(2 * i, 0, a.get(0, 0));
  56.                         m.set(2 * i + 1, 0, a.get(1, 0));
  57.                         i++;
  58.                 }
  59.                 return m;
  60.  
  61.         }
  62.  
  63.         private ArrayList<Matrix> matrixToArray(Matrix matrix) {
  64.                 ArrayList<Matrix> array = new ArrayList<Matrix>();
  65.  
  66.                 for (int i = 0; i < matrix.getRowDimension(); i++) {
  67.                         Matrix m = new Matrix(2, 1);
  68.                         m.set(0, 0, matrix.get(i, 0));
  69.                         m.set(1, 0, matrix.get(++i, 0));
  70.                         array.add(m);
  71.                 }
  72.  
  73.                 return array;
  74.  
  75.         }
  76.  
  77.         public ArrayList<Matrix> startInteriorPoint(int alghoritm) {
  78.                 interiorAll = 0;
  79.                 newtonAv = 0.;
  80.                 newtonAll = 0;
  81.  
  82.                 int n = buildings.size();
  83.                 Matrix temp = arrayToMatrix(buildings);
  84.                 double t = t0;
  85.                 while (((double) n) / t > epsilon) {
  86.                         interiorAll++;
  87.                         t = t * mi;
  88.                         temp = newton(temp, t, alghoritm);
  89.                 }
  90.                 buildings = matrixToArray(temp);
  91.                 // for(Matrix b: buildings){
  92.                 // System.out.println(b.get(0,0)+" "+b.get(1,0));
  93.                 // }
  94.                 newtonAv = (double) newtonAll / interiorAll;
  95.                 interiorResult = o(temp);
  96.                 return buildings;
  97.         }
  98.  
  99.         public double o(Matrix arg) {
  100.                 int p = 2;
  101.                 double res = 0D;
  102.                 for (int i = 0; i < this.supplies.size(); i++) {
  103.                         Zasob o = this.supplies.get(i);
  104.                         double norm = Math.sqrt(normSquare(arg.get(p, 0),
  105.                                         arg.get(p + 1, 0), o.getX(), o.getY()));
  106.                         res += o.getA() * g(norm - o.getR());
  107.                         p += 2;
  108.                 }
  109.                 return res;
  110.         }
  111.  
  112.         private double g(double x) {
  113.                 return 0.5D * (x + Math.sqrt(x * x + gamma));
  114.         }
  115.  
  116.         private double minimize(int alg, Matrix x, Matrix dX, double t) {
  117.                 switch (alg) {
  118.                 case BACKTRACKING:
  119.                         return backtracking(x, dX, t);
  120.                 case GOLDENRATIO:
  121.                         return goldenRatio(x, dX, t);
  122.                 default:
  123.                         return 0;
  124.                 }
  125.         }
  126.  
  127.         public Matrix newton(Matrix start, double t, int alghoritm) {
  128.                 Matrix temp = (Matrix) start.clone();
  129.                 double gradNorm;
  130.                 double ni;
  131.                 do {
  132.                         Matrix gradient = gradient(temp, t);
  133.                         Matrix hessian = hessian(temp, t);
  134.                         Matrix hessianInv = hessian.inverse();
  135.                         Matrix hessianInvXgradient = hessianInv.times(gradient);
  136.                         gradNorm = gradient.transpose().times(hessianInvXgradient)
  137.                                         .get(0, 0);
  138.                         hessianInvXgradient = hessianInvXgradient.times(-1.0);
  139.  
  140.                         ni = minimize(alghoritm, temp, hessianInvXgradient, t);
  141.                         temp = temp.plus(hessianInvXgradient.times(ni));
  142.  
  143.                 } while (gradNorm >= epsilon && ni >= 1.0E-8D && newtonAll++ < 10000);
  144.  
  145.                 return temp;
  146.  
  147.         }
  148.  
  149.         private double gt(Matrix x, double t) {
  150.  
  151.                 double ogran = 0D;
  152.                 int p = 2;
  153.                 for (int i = 0; i < supplies.size(); i++) {
  154.                         double normKwa = normSquare(x.get(p, 0), x.get(p + 1, 0),
  155.                                         x.get(0, 0), x.get(1, 0));
  156.                         ogran += Math.log(r * r - normKwa);
  157.                         p += 2;
  158.                 }
  159.                 return t * o(x) - ogran;
  160.  
  161.         }
  162.  
  163.         private double normSquare(double x1, double x2, double y1, double y2) {
  164.                 return (x1 - y1) * (x1 - y1) + (x2 - y2) * (x2 - y2);
  165.         }
  166.  
  167.         private double backtracking(Matrix start, Matrix dX, double t) {
  168.                 Matrix gradient = gradient(start, t).transpose();
  169.                 Matrix step = gradient.times(dX);
  170.                 double ni = 1;
  171.                 double newVal;
  172.                 do {
  173.                         ni *= beta;
  174.                         Matrix newStart = start.plus(dX.times(ni));
  175.                         newVal = gt(newStart, t);
  176.                 } while ((Double.isNaN(newVal))
  177.                                 || (newVal > gt(start, t) + alfa * ni * step.get(0, 0)));
  178.                 System.out.println(ni);
  179.                 return ni;
  180.         }
  181.  
  182.         private double goldenRatio(Matrix start, Matrix dX, double t) {
  183.                 double alfa = 0.618;
  184.                 double a = 0;
  185.                 double b = 1;
  186.  
  187.                 double norm = dX.norm2();
  188.  
  189.                 double lambda = a + (1 - alfa) * (b - a);
  190.                 double mi = a + alfa * (b - a);
  191.  
  192.                 while ((b - a) * norm > epsilon) {
  193.                         double left = gt(start.plus(dX.times(lambda)), t);
  194.                         double right = gt(start.plus(dX.times(mi)), t);
  195.                         if (left > right) {
  196.                                 a = lambda;
  197.                                 lambda = mi;
  198.                                 mi = a + alfa * (b - a);
  199.                         } else {
  200.                                 b = mi;
  201.                                 mi = lambda;
  202.                                 lambda = a + (1 - alfa) * (b - a);
  203.                         }
  204.                 }
  205.  
  206.                 return a / 2 + alfa * b / 2;
  207.         }
  208.  
  209.         private Matrix gradient(Matrix m, double t) {
  210.                 int n = buildings.size();
  211.                 Matrix res = new Matrix(2 * n, 1);
  212.  
  213.                 double zeroX = 0D;
  214.                 double zeroY = 0D;
  215.                 int p = 2;
  216.                 for (int i = 0; i < n - 1; i++) {
  217.                         double norm2 = normSquare(m.get(p, 0), m.get(p + 1, 0),
  218.                                         m.get(0, 0), m.get(1, 0));
  219.                         zeroX += 2.0D * (m.get(p, 0) - m.get(0, 0)) / (r * r - norm2);
  220.                         zeroY += 2.0D * (m.get(p + 1, 0) - m.get(1, 0)) / (r * r - norm2);
  221.                         p += 2;
  222.                 }
  223.                 res.set(0, 0, -zeroX);
  224.                 res.set(1, 0, -zeroY);
  225.  
  226.                 p = 2;
  227.                 for (int i = 0; i < n - 1; i++) {
  228.                         Zasob o = supplies.get(i);
  229.                         double norm2X = normSquare(m.get(p, 0), m.get(p + 1, 0),
  230.                                         m.get(0, 0), m.get(1, 0));
  231.                         double norm2Z = normSquare(m.get(p, 0), m.get(p + 1, 0), o.getX(),
  232.                                         o.getY());
  233.                         double sqrtnorm2Z = Math.sqrt(norm2Z);
  234.  
  235.                         double x1 = t * o.getA() * g1(sqrtnorm2Z - o.getR()) / sqrtnorm2Z;
  236.  
  237.                         double x2 = 2.0D / (r * r - norm2X);
  238.  
  239.                         double X = x1 * (m.get(p, 0) - o.getX()) + x2
  240.                                         * (m.get(p, 0) - m.get(0, 0));
  241.                         double Y = x1 * (m.get(p + 1, 0) - o.getY()) + x2
  242.                                         * (m.get(p + 1, 0) - m.get(1, 0));
  243.  
  244.                         res.set(p, 0, X);
  245.                         res.set(p + 1, 0, Y);
  246.  
  247.                         p += 2;
  248.                 }
  249.  
  250.                 return res;
  251.  
  252.         }
  253.  
  254.         private double g1(double x) {
  255.                 return 0.5 * (x / Math.sqrt(x * x + gamma)) + 0.5;
  256.         }
  257.  
  258.         public Matrix hessian(Matrix arg, double t) {
  259.                 Matrix res = new Matrix(2 * this.supplies.size() + 2,
  260.                                 2 * this.supplies.size() + 2);
  261.  
  262.                 double zeroXX = 0D;
  263.                 double zeroYY = 0D;
  264.                 double zeroXY = 0D;
  265.                 int p = 2;
  266.                 for (int i = 0; i < this.supplies.size(); i++) {
  267.                         double normKwa = r
  268.                                         * r
  269.                                         - normSquare(arg.get(p, 0), arg.get(p + 1, 0),
  270.                                                         arg.get(0, 0), arg.get(1, 0));
  271.                         double currXX = 4.0D / (normKwa * normKwa)
  272.                                         * (arg.get(p, 0) - arg.get(0, 0))
  273.                                         * (arg.get(p, 0) - arg.get(0, 0)) + 2.0D / normKwa;
  274.                         double currYY = 4.0D / (normKwa * normKwa)
  275.                                         * (arg.get(p + 1, 0) - arg.get(1, 0))
  276.                                         * (arg.get(p + 1, 0) - arg.get(1, 0)) + 2.0D / normKwa;
  277.                         double currXY = 4.0D / (normKwa * normKwa)
  278.                                         * (arg.get(p, 0) - arg.get(0, 0))
  279.                                         * (arg.get(p + 1, 0) - arg.get(1, 0));
  280.  
  281.                         res.set(p, 0, -currXX);
  282.                         res.set(0, p, -currXX);
  283.  
  284.                         res.set(p + 1, 1, -currYY);
  285.                         res.set(1, p + 1, -currYY);
  286.  
  287.                         res.set(p, 1, -currXY);
  288.                         res.set(1, p, -currXY);
  289.                         res.set(p + 1, 0, -currXY);
  290.                         res.set(0, p + 1, -currXY);
  291.  
  292.                         zeroXX += currXX;
  293.                         zeroYY += currYY;
  294.                         zeroXY += currXY;
  295.  
  296.                         p += 2;
  297.                 }
  298.                 res.set(0, 0, zeroXX);
  299.                 res.set(1, 1, zeroYY);
  300.                 res.set(1, 0, zeroXY);
  301.                 res.set(0, 1, zeroXY);
  302.  
  303.                 p = 2;
  304.  
  305.                 for (int i = 0; i < this.supplies.size(); i++) {
  306.                         Zasob o = (Zasob) this.supplies.get(i);
  307.                         double norm2X = normSquare(arg.get(p, 0), arg.get(p + 1, 0),
  308.                                         arg.get(0, 0), arg.get(1, 0));
  309.                         double norm2Z = normSquare(arg.get(p, 0), arg.get(p + 1, 0),
  310.                                         o.getX(), o.getY());
  311.                         double sqrtnorm2Z = Math.sqrt(norm2Z);
  312.                         double cubenormZ = sqrtnorm2Z * sqrtnorm2Z * sqrtnorm2Z;
  313.  
  314.                         double x1 = t
  315.                                         * o.getA()
  316.                                         * (g2(sqrtnorm2Z - o.getR()) * sqrtnorm2Z - g1(sqrtnorm2Z
  317.                                                         - o.getR())) / cubenormZ;
  318.                         double x2 = t * o.getA() * g1(sqrtnorm2Z - o.getR()) / sqrtnorm2Z;
  319.                         double x3 = 4.0D / ((r * r - norm2X) * (r * r - norm2X));
  320.                         double x4 = 2.0D / (r * r - norm2X);
  321.  
  322.                         double wartXX = x1 * (arg.get(p, 0) - o.getX())
  323.                                         * (arg.get(p, 0) - o.getX()) + x2 + x3
  324.                                         * (arg.get(p, 0) - arg.get(0, 0))
  325.                                         * (arg.get(p, 0) - arg.get(0, 0)) + x4;
  326.                         double wartYY = x1 * (arg.get(p + 1, 0) - o.getY())
  327.                                         * (arg.get(p + 1, 0) - o.getY()) + x2 + x3
  328.                                         * (arg.get(p + 1, 0) - arg.get(1, 0))
  329.                                         * (arg.get(p + 1, 0) - arg.get(1, 0)) + x4;
  330.                         double wartXY = x1 * (arg.get(p, 0) - o.getX())
  331.                                         * (arg.get(p + 1, 0) - o.getY()) + x3
  332.                                         * (arg.get(p, 0) - arg.get(0, 0))
  333.                                         * (arg.get(p + 1, 0) - arg.get(1, 0));
  334.  
  335.                         res.set(p, p, wartXX);
  336.                         res.set(p + 1, p + 1, wartYY);
  337.                         res.set(p, p + 1, wartXY);
  338.                         res.set(p + 1, p, wartXY);
  339.  
  340.                         p += 2;
  341.                 }
  342.  
  343.                 return res;
  344.         }
  345.  
  346.         private double g2(double x) {
  347.                 double sqrt = Math.sqrt(x * x + gamma);
  348.                 return 0.5 * (1.0 / sqrt - x * x / (sqrt * sqrt * sqrt));
  349.         }
  350. }
create a new version of this paste RAW Paste Data