Andrey Akinshin Add prony-2d-geometry  over 8 years ago

Commit id: d602fbddfc9a4552be7490c57f198373f4fb392d

deletions | additions      

         

\subsection{Real Prony, Geometry (WIP)}        Binary files /dev/null and b/figures/prony-2d-geometry/plot2d-h0.025.png differ       Binary files /dev/null and b/figures/prony-2d-geometry/plot2d-h0.05.png differ       Binary files /dev/null and b/figures/prony-2d-geometry/plot2d-h0.1.png differ       Binary files /dev/null and b/figures/prony-2d-geometry/plot2d-h0.2.png differ       Binary files /dev/null and b/figures/prony-2d-geometry/plot2d-h0.4.png differ       Binary files /dev/null and b/figures/prony-2d-geometry/plot3d-h0.025.png differ       Binary files /dev/null and b/figures/prony-2d-geometry/plot3d-h0.05.png differ       Binary files /dev/null and b/figures/prony-2d-geometry/plot3d-h0.1.png differ       Binary files /dev/null and b/figures/prony-2d-geometry/plot3d-h0.2.png differ       Binary files /dev/null and b/figures/prony-2d-geometry/plot3d-h0.4.png differ          

height = 700\nwidth = 500         

figures/prony-2d-countour3d/prony-2d-countour3d.gif  Numerical_Prony2D_TrigonometricReal2.tex  Numerical_Prony2D_TrigonometricRealNoisy.tex  figures/prony-2d-noisy/prony-2d-noisy.png Numerical_Prony2D_Geometry.tex  figures/prony-2d-geometry/plot2d-h0.025.png  figures/prony-2d-geometry/plot3d-h0.025.png  figures/prony-2d-geometry/plot2d-h0.05.png  figures/prony-2d-geometry/plot3d-h0.05.png  figures/prony-2d-geometry/plot2d-h0.1.png  figures/prony-2d-geometry/plot3d-h0.1.png  figures/prony-2d-geometry/plot2d-h0.2.png  figures/prony-2d-geometry/plot3d-h0.2.png  figures/prony-2d-geometry/plot2d-h0.4.png  figures/prony-2d-geometry/plot3d-h0.4.png           

library(rgl)  library(ggplot2)  library(dplyr)  # generator  system("mcs -o+ prony-2d-geometry.cs")  generate <- function(h, make.2d = T, make.3d = T, show = T, save = F) {  # create data  system(paste("prony-2d-geometry.exe", h))    # loading and processing  df <- read.csv("data.csv", sep=";")  breaks <- c(0, 0.0001, 0.0003, 0.0005, 0.0007, 0.001, 0.002, 0.006, 0.01, 0.02, 0.06, 100)  clen <- length(levels(cut(df$value, breaks = breaks)))  colorlut <- rainbow(clen)  col <- colorlut[df$category]  df <- df %>%   mutate(category = as.numeric(cut(value, breaks = breaks, labels = 1:clen))) %>%   mutate(interval = cut(value, breaks = breaks)) %>%   mutate(color = colorlut[category])    # plots  if (show | save) {  if (make.2d) {  p <- ggplot(na.omit(df), aes(x = x, y = y, fill = interval)) + geom_raster() +  ggtitle(paste("h =", h)) + xlab("x") + ylab("y") + guides(fill = guide_legend(title = "Eps")) +  scale_fill_manual(values = colorlut) + theme_bw()  if (show)  print(p)  if (save)  ggsave(paste0("../figures/prony-2d-geometry/plot2d-h", h, ".png"), p)  }  if (make.3d) {  xs <- unique(df$x)  ys <- unique(df$y)  open3d(windowRect = c(-1000, 100, -500, 600))  surface3d(xs, ys, matrix(df$value, nrow=length(xs), ncol=length(ys)), color = df$color, back = "lines")  par3d(userMatrix = rotationMatrix(-pi/4 - pi, -0.15, 0.3, 1))  if (save)  snapshot3d(paste0("../figures/prony-2d-geometry/plot3d-h", h, ".png"))  if (!show)  rgl.close()  }  }  }  for (h in c(0.025, 0.05, 0.1, 0.2, 0.4))  generate(h, show = F, save = T)  # generate(0.1)           

using System;  using System.Diagnostics;  using System.Globalization;  using System.IO;  using System.Linq;  using System.Runtime.CompilerServices;  using System.Text;  public class Program  {  public static void Main(string[] args)  {  var h = args.Length == 0 ? 0.1 : double.Parse(args[0], CultureInfo.InvariantCulture);  var sw = Stopwatch.StartNew();  var builder = new StringBuilder();  builder.AppendLine("x;y;value");  Calc(builder, h);  File.WriteAllText("data.csv", builder.ToString());  sw.Stop();  Console.WriteLine($"Done, TotalTime = {sw.ElapsedMilliseconds} ms");  }  private static void Calc(StringBuilder builder, double h)  {  double[] xs = Seq(-h, h, 501), ys = Seq(0, 2 * h, 501), das = Seq(-0.01, 0.01, 51);  double refm1 = h, refm2 = h * h, refm3 = h * h * h;  foreach (var x in xs)  foreach (var y in ys)  {  var best = 1e100;  var a0 = 2 - Clamp(Math.Abs(x - y) < 1e-9 ? 0 : (h - 2 * x) / (y - x), 0, 2);  foreach (var da in das)  {  var a = Clamp(a0 + da, 0, 2);  var b = 2 - a;  var m1 = a * x + b * y;  var m2 = a * x * x + b * y * y;  var m3 = a * x * x * x + b * y * y * y;  var d = Sqr(refm1 - m1) + Sqr(refm2 - m2) + Sqr(refm3 - m3);  best = Math.Min(d, best);  }  builder.AppendLine($"{x:0.0000};{y:0.0000};{Math.Sqrt(best):0.0000}");  }  }  private static double[] Seq(double from, double to, double length) =>  Enumerable.Range(0, (int)length).Select(i => @from + (to - @from) / ((int)length - 1) * i).ToArray();  [MethodImpl(MethodImplOptions.AggressiveInlining)]  public static double Sqr(double value) => value * value;  [MethodImpl(MethodImplOptions.AggressiveInlining)]  public static double Clamp(double value, double min, double max) => Math.Min(max, Math.Max(min, value));  }