Merge branch 'master' of https://github.com/grasingerm/polymer-stats

2D/inc/acceptance.jl

@@ -0,0 +1,39 @@
+@inline function metropolis_acc(kT::Real, dU::Real, ϵ::Real)
+  return ϵ <= (exp(-dU / kT));
+end
+
+function kawasaki_acc(kT::Real, dU::Real, ϵ::Real)
+  boltz = exp(-dU / (2*kT));
+  anti_boltz = exp(dU / (2*kT));
+  return ( ϵ <= (boltz / (boltz + anti_boltz)) );
+end
+
+abstract type Acceptor end;
+
+mutable struct Metropolis <: Acceptor
+  logπ_prev::Float64;
+  weight_function::WeightFunction;
+end
+
+logπ_chain(chain::EAPChain) = -chain.U / chain.kT;
+
+function logπ_chain(chain::EAPChain, log_wf::WeightFunction)
+  return -chain.U / chain.kT + log_wf(chain);
+end
+
+function Metropolis(chain::EAPChain, wf::WeightFunction)
+  return Metropolis(logπ_chain(chain, wf), wf);
+end
+
+# alpha is a ratio of trial move probabilities (related to clustering)
+function (metro::Metropolis)(chain::EAPChain, ϵ::Real; α::Real = 1.0)
+  metro.weight_function(chain);
+  logπ_curr = logπ_chain(chain, metro.weight_function) + log(α);
+  if (logπ_curr >= metro.logπ_prev) || (ϵ < exp(logπ_curr - metro.logπ_prev))
+    metro.logπ_prev = logπ_curr;
+    return true;
+  else
+    return false;
+  end
+  return false;
+end

2D/inc/algorithm_help_msg.jl

@@ -0,0 +1,28 @@
+ _ALGO_HELP_MSG_ = """algorithm (ISRES | PRAXIS | SBPLX | SLSQP | LBFGS | AUGLAG | Newton | Anneal)
+  
+  Separate algorithms by commas to run them in succession; e.g. LBGFS,Newton
+
+  Algorithms:
+  ISRES   -   Improved stochastic ranking evolution strategy; 
+              nonlinear rained global optimization;
+              note: requires arbitrary bounds on unknowns
+  PRAXIS  -   PRincipal AXIS; gradient-free local optimization
+  SBPLX   -   Subplex; variant of Nelder-Mead that is more efficient and robust;
+              gradient-free local optimization
+  SLSQP   -   Sequential quadratic programming for nonlineraly rained,
+              gradient-based optimizationa
+  LBFGS   -   low-storage BFGS; gradient-based unrained optimization
+  AUGLAG  -   Augmented Lagrangian; rained optimization with LBFGS
+  Newton  -   My implementation of the Newton-Raphson method
+  Anneal  -   My implementation of simulated annealing
+
+  Special syntax: try algorithm1,algorithm2 then backup_algorithms
+  This first tries an algorithm that may be less robust but faster and
+  if it fails then an algorithm that is more robust and slower.
+
+  Note: in general, I've found SBPLX to be incredibly robust (for this problem).
+  If you have no clue at an initial guess, my advice is to start with some
+  simulated anneeling, hit with the SBPLX, and then finish it off with a little
+  of Newton's method. I've tested this with intentionally pathological initial
+  guesses and often converge to the solution anyway :)
+""";

2D/inc/average.jl

@@ -0,0 +1,124 @@
+
+
+const Accessor = Function;
+
+abstract type Averager end;
+
+# TODO: consider creating an averager which uses arbitrary precision arithmetic
+mutable struct StandardAverager{T,N} <: Averager
+  value::T;
+  normalizer::N;
+  accessor::Accessor;
+end
+
+function StandardAverager(accessor::Accessor, chain::EAPChain)
+  return StandardAverager(
+                          0.0 * accessor(chain),
+                          0.0,
+                          accessor
+                         );
+end
+
+function StandardAverager(dt::DataType, accessor::Accessor, chain::EAPChain)
+  return StandardAverager(
+                          dt(0.0 * accessor(chain)),
+                          dt(0.0),
+                          accessor
+                         );
+end
+
+function StandardAverager(dt1::DataType, dt2::DataType, accessor::Accessor, chain::EAPChain)
+  return StandardAverager(
+                          dt2(0.0 * accessor(chain)),
+                          dt1(0.0),
+                          accessor
+                         );
+end
+
+@inline get_avg(stdavg::StandardAverager) = stdavg.value / stdavg.normalizer;
+
+function record!(stdavg::StandardAverager, chain::EAPChain)
+  stdavg.value += stdavg.accessor(chain);
+  stdavg.normalizer += 1;
+end
+
+function record!(stdavg::StandardAverager{<:Vector}, chain::EAPChain)
+  stdavg.value[:] += stdavg.accessor(chain);
+  stdavg.normalizer += 1;
+end
+
+abstract type WeightFunction end;
+
+struct UmbrellaAverager{T,N} <: Averager
+  stdavg::StandardAverager{T,N};
+  weightFunction::WeightFunction;
+end
+
+function UmbrellaAverager(ac::Accessor, wf::WeightFunction, chain::EAPChain)
+  return UmbrellaAverager(StandardAverager(ac, chain), wf);
+end
+
+@inline get_avg(ua::UmbrellaAverager) = get_avg(ua.stdavg);
+
+function record!(ua::UmbrellaAverager, chain::EAPChain)
+  expw = exp(ua.weightFunction(chain));
+  ua.stdavg.value += ua.stdavg.accessor(chain) / expw;
+  ua.stdavg.normalizer += 1.0 / expw;
+end
+
+function record!(ua::UmbrellaAverager{Vector{<:Real}}, chain::EAPChain)
+  expw = exp(ua.weightFunction(chain));
+  ua.stdavg.value[:] += ua.stdavg.accessor(chain) / expw;
+  ua.stdavg.normalizer += 1.0 / expw;
+end
+
+function record!(ua::UmbrellaAverager{Dec128}, chain::EAPChain)
+  expw = exp(Dec128(ua.weightFunction(chain)));
+  ua.stdavg.value += ua.stdavg.accessor(chain) / expw;
+  ua.stdavg.normalizer += 1.0 / expw;
+end
+
+function record!(ua::UmbrellaAverager{Vector{Dec128}}, chain::EAPChain)
+  expw = exp(Dec128(ua.weightFunction(chain)));
+  ua.stdavg.value[:] += ua.stdavg.accessor(chain) / expw;
+  ua.stdavg.normalizer += 1.0 / expw;
+end
+
+function record!(ua::UmbrellaAverager{BigFloat}, chain::EAPChain)
+  expw = exp(BigFloat(ua.weightFunction(chain)));
+  ua.stdavg.value += ua.stdavg.accessor(chain) / expw;
+  ua.stdavg.normalizer += 1.0 / expw;
+end
+
+function record!(ua::UmbrellaAverager{Vector{BigFloat}}, chain::EAPChain)
+  expw = exp(BigFloat(ua.weightFunction(chain)));
+  ua.stdavg.value[:] += ua.stdavg.accessor(chain) / expw;
+  ua.stdavg.normalizer += 1.0 / expw;
+end
+
+struct WeightlessFunction <: WeightFunction
+end
+
+(wf::WeightlessFunction)(::EAPChain) = 1.0;
+
+struct AntiDipoleWeightFunction <: WeightFunction
+  log_gauge::Float64;
+end
+
+# forward call
+@inline AntiDipoleWeightFunction(chain::EAPChain) = AntiDipoleWeightFunction(chain, chain.μ);
+function AntiDipoleWeightFunction(chain::EAPChain, dr::DielectricResponse)
+  return AntiDipoleWeightFunction(-(dr.K1 + 2 * dr.K2) * chain.E0*chain.E0 *
+                                   n(chain) / (3 * chain.kT));
+end
+
+function AntiDipoleWeightFunction(chain::EAPChain, dr::PolarResponse)
+  return AntiDipoleWeightFunction(-eigvals(dr.M)[end] * chain.E0 *
+                                   n(chain) / (3 * chain.kT));
+end
+
+function (wf::AntiDipoleWeightFunction)(chain::EAPChain)
+  return (sum(chain.us) / chain.kT * 
+          (0.2 + 0.8*(exp(-(chain.Fx^2+chain.Fz^2)/chain.kT))) 
+          - wf.log_gauge);
+end

2D/inc/dielectric_dipole.jl

@@ -0,0 +1,4 @@
+function μ(E0::Real, K1::Real, K2::Real, cϕ::Real, sϕ::Real)
+  n̂i = n̂(cϕ, sϕ);
+  return ((K1 - K2) * E0 * cϕ * n̂i) + (K2 * [E0; 0.0]);
+end

2D/inc/dipole_response.jl

@@ -0,0 +1,27 @@
+include(joinpath(@__DIR__, "types.jl"));
+
+abstract type DipoleResponse end
+
+(dr::DipoleResponse)(E0::Real) = begin; error("Not yet implemented"); E0; end
+
+function μ_dielectric(E0::Real, K1::Real, K2::Real, cϕ::Real, sϕ::Real)
+  n̂i = n̂(cϕ, sϕ);
+  return ((K1 - K2) * E0 * cϕ * n̂i) + (K2 * [E0; 0.0]);
+end
+
+struct DielectricResponse <: DipoleResponse
+  K1::Real;
+  K2::Real;
+end
+
+function (dr::DielectricResponse)(E0::Real, cϕ::Real, sϕ::Real)
+  return μ_dielectric(E0, dr.K1, dr.K2, cϕ, sϕ);
+end
+
+struct PolarResponse <: DipoleResponse
+  M::FdMatrix;
+end
+
+function (pr::PolarResponse)(::Real, cϕ::Real, sϕ::Real)
+  return pr.M * n̂(cϕ, sϕ);
+end