From 1932b4f0b0d9ce8ec53cdee0fc70b6feba67a845 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Fri, 24 Sep 2021 13:55:26 -0400 Subject: [PATCH 01/17] tested clustering algorithm against Langevin force-stretch relationship --- inc/acceptance.jl | 5 +- inc/eap_chain.jl | 107 +++++++++++- inc/energy.jl | 7 + mcmc_clustering_eap_chain.jl | 328 +++++++++++++++++++++++++++++++++++ mcmc_eap_chain.jl | 10 +- 5 files changed, 452 insertions(+), 5 deletions(-) create mode 100644 mcmc_clustering_eap_chain.jl diff --git a/inc/acceptance.jl b/inc/acceptance.jl index 1b36f40..c0d0d2a 100644 --- a/inc/acceptance.jl +++ b/inc/acceptance.jl @@ -25,9 +25,10 @@ function Metropolis(chain::EAPChain, wf::WeightFunction) return Metropolis(logπ_chain(chain, wf), wf); end -function (metro::Metropolis)(chain::EAPChain, ϵ::Real) +# 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π_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; diff --git a/inc/eap_chain.jl b/inc/eap_chain.jl index 43b83ee..4320b8e 100644 --- a/inc/eap_chain.jl +++ b/inc/eap_chain.jl @@ -8,6 +8,7 @@ include(joinpath(@__DIR__, "dipole_response.jl")); θ_dist = Uniform(0.0, π); # the way that the energy functions are organized is a mess; I'm sorry +# also... not sorry abstract type Energy end mutable struct EAPChain @@ -89,6 +90,8 @@ function EAPChain(pargs::Dict) NonInteractingEnergy(); elseif pargs["energy-type"] == "interacting" InteractingEnergy(); + elseif pargs["energy-type"] == "Ising" + IsingEnergy(); else error("energy-type is not understood."); end, @@ -158,13 +161,30 @@ function U_interaction(chain::EAPChain) r̂ = r / rmag; r3 = r2*rmag; μi = view(chain.μs, :, i); - μj = view(chain.μs, :, i); - U += (dot(μi, μj) - 3*dot(μi, r̂)*dot(μj, r̂)) / (4*π); + μj = view(chain.μs, :, j); + U += (dot(μi, μj) - 3*dot(μi, r̂)*dot(μj, r̂)) / (4*π*r3); end end return U; end +# TODO: consider rewriting this to make better use of past calculations +# this is probably by far the slowest calculation *shrug emoji* +function U_Ising(chain::EAPChain) + U = 0.0; + @inbounds for i=1:n(chain)-1 + r = chain.xs[:, i] - chain.xs[:, i+1]; + r2 = dot(r, r); + rmag = sqrt(r2); + r̂ = r / rmag; + r3 = r2*rmag; + μi = view(chain.μs, :, i); + μj = view(chain.μs, :, i+1); + U += (dot(μi, μj) - 3*dot(μi, r̂)*dot(μj, r̂)) / (4*π*r3); + end + return U; +end + function move!(chain::EAPChain, idx::Int, dϕ::Real, dθ::Real) @inbounds begin chain.ϕs[idx] += dϕ; @@ -189,6 +209,89 @@ function move!(chain::EAPChain, idx::Int, dϕ::Real, dθ::Real) return true; end +function flip_n!(chain::EAPChain, idx::Int) + move!(chain, idx, π, π - 2*chain.θs[idx]) +end + +pflip_linear(ip::Real) = (1 + ip) / 2; + +function cluster_flip!(chain::EAPChain, idx::Int; + pflip::Function = pflip_linear, + ϵflip::Real = 1.0, + ηerr::Real = 1e-10 + ) + + # grow the cluster to the right + upper_p = NaN; + upper_idx = idx; + while true + if upper_idx >= n(chain) + upper_p = 0.0; + break; + end + upper_p = pflip_linear(dot(n̂j(chain, upper_idx), n̂j(chain, upper_idx+1))); + if rand() <= upper_p + upper_idx += 1 + else + break; + end + end + @assert(!isnan(upper_p) && upper_idx <= n(chain)); + + # grow the cluster to the left + lower_p = NaN; + lower_idx = idx; + while true + if lower_idx <= 1 + lower_p = 0.0; + break; + end + lower_p = pflip_linear(dot(n̂j(chain, lower_idx), n̂j(chain, lower_idx-1))); + if rand() <= lower_p + lower_idx -= 1 + else + break; + end + end + @assert(!isnan(lower_p) && lower_idx >= 1); + + if rand() <= ϵflip + prev_n̂s = copy(chain.n̂s[:, lower_idx:upper_idx]); + for i in lower_idx:upper_idx + flip_n!(chain, i); + end + + if norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) > ηerr + println("prev"); + display(prev_n̂s) + println(); + println("curr"); + display(chain.n̂s[:, lower_idx:upper_idx]) + println(); + @show(norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf)) + end + @assert(norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) < ηerr); + + new_upper_p = if upper_idx < n(chain) + pflip_linear(dot(n̂j(chain, upper_idx), n̂j(chain, upper_idx+1))); + else + 0 + end + new_lower_p = if lower_idx > 1 + pflip_linear(dot(n̂j(chain, lower_idx), n̂j(chain, lower_idx-1))); + else + 0 + end + return ( + ((1 - new_upper_p)*(1 - new_lower_p)) / + ((1 - upper_p)*(1 - lower_p)) + ); + else + return 1.0; + end + +end + function move!(chain::EAPChain, idx::Int, dϕ::Real, dθ::Real, r0::AbstractVector; frac_mv::Real = 0.15, acc_tol::Real = 1e-1) diff --git a/inc/energy.jl b/inc/energy.jl index b0b2cd3..b055910 100644 --- a/inc/energy.jl +++ b/inc/energy.jl @@ -14,3 +14,10 @@ struct InteractingEnergy <: Energy end return (sum(chain.us) + U_interaction(chain) - dot(end_to_end(chain), [chain.Fx; 0.0; chain.Fz])); end + +struct IsingEnergy <: Energy end + +@inline function (::IsingEnergy)(chain::EAPChain) + return (sum(chain.us) + U_Ising(chain) + - dot(end_to_end(chain), [chain.Fx; 0.0; chain.Fz])); +end diff --git a/mcmc_clustering_eap_chain.jl b/mcmc_clustering_eap_chain.jl new file mode 100644 index 0000000..430f0fb --- /dev/null +++ b/mcmc_clustering_eap_chain.jl @@ -0,0 +1,328 @@ +using ArgParse; +using Distributions; +using LinearAlgebra; +using Logging; +using DelimitedFiles; +using ProfileView; +using DecFP; +using Quadmath; + +include(joinpath("inc", "eap_chain.jl")); +include(joinpath("inc", "average.jl")); +include(joinpath("inc", "acceptance.jl")); + +s = ArgParseSettings(); +@add_arg_table! s begin + "--E0", "-e" + help = "magnitude of electric field" + arg_type = Float64 + default = 0.0 + "--chain-type", "-T" + help = "chain type (dielectric|polar)" + arg_type = String + default = "dielectric" + "--K1", "-J" + help = "dipole susceptibility along the monomer axis (dielectric chain)" + arg_type = Float64 + default = 1.0 + "--K2", "-K" + help = "dipole susceptibility orthogonal to the monomer axis (dielectric chain)" + arg_type = Float64 + default = 0.0 + "--mu", "-m" + arg_type = Float64 + default = 1e-2 + help = "dipole magnitude (electret chain)" + "--energy-type", "-u" + help = "energy type (noninteracting|interacting|Ising)" + arg_type = String + default = "noninteracting" + "--kT", "-k" + help = "dimensionless temperature" + arg_type = Float64 + default = 1.0 + "--Fz", "-F" + help = "force in the z-direction (direction of E-field; force ensemble)" + arg_type = Float64 + default = 0.0 + "--Fx", "-G" + help = "force in the x-direction (force ensemble)" + arg_type = Float64 + default = 0.0 + "--mlen", "-b" + help = "monomer length" + arg_type = Float64 + default = 1.0 + "--num-monomers", "-n" + help = "number of monomers" + arg_type = Int + default = 100 + "--num-steps", "-N" + help = "number of steps" + arg_type = Int + default = convert(Int, 1e6) + "--phi-step", "-p" + help = "maximum ϕ step length" + arg_type = Float64; + default = 3*π / 8; + "--theta-step", "-q" + help = "maximum θ step length" + arg_type = Float64; + default = 3*π / 16; +# "--do-clustering" +# help = "trial moves with clustering" +# action = :store_true + "--cluster-prob" + help = "probability of flipping a cluster" + arg_type = Float64; + default = 0.5; + "--step-adjust-lb", "-L" + help = "adjust step sizes if acc. ratio below this threshold" + arg_type = Float64 + default = 0.15 + "--step-adjust-ub", "-U" + help = "adjust step sizes if acc. ratio above this threshold" + arg_type = Float64 + default = 0.40 + "--step-adjust-scale", "-A" + help = "scale factor for adjusting step sizes (> 1.0)" + arg_type = Float64 + default = 1.1 + "--steps-per-adjust", "-S" + help = "steps between step size adjustments" + arg_type = Int + default = 2500 + "--umbrella-sampling", "-B" + help = "use umbrella sampling (w/ electrostatic weight function)" + action = :store_true + "--update-freq" + help = "update frequency (seconds)" + arg_type = Float64; + default = 15.0; + "--verbose", "-v" + help = "verbosity level: 0-nothing, 1-errors, 2-warnings, 3-info" + arg_type = Int + default = 3 + "--prefix", "-P" + help = "prefix for output files" + arg_type = String + default = "eap-mcmc" + "--postfix", "-Q" + help = "postfix for output files" + arg_type = String + default = "" + "--stepout", "-s" + help = "steps between storing microstates" + arg_type = Int + default = 500 + "--numeric-type" + help = "numerical data type for averaging (float64|float128|dec128|big)" + arg_type = String + default = "float64" + "--profile", "-Z" + help = "profile the program" + action = :store_true +end + +pargs = parse_args(s); + +if pargs["verbose"] == 3 + global_logger(ConsoleLogger(stderr, Logging.Info)); +elseif pargs["verbose"] == 2 + global_logger(ConsoleLogger(stderr, Logging.Warn)); +elseif pargs["verbose"] == 1 + global_logger(ConsoleLogger(stderr, Logging.Error)); +else + global_logger(Logging.NullLogger()); +end + +function mcmc(nsteps::Int, pargs) + ϕstep, θstep = pargs["phi-step"], pargs["theta-step"]; + dϕ_dist = Uniform(-ϕstep, ϕstep); + dθ_dist = Uniform(-θstep, θstep); + chain = EAPChain(pargs); + chain.U = U(chain); + wf = AntiDipoleWeightFunction(chain); + acceptor = Metropolis( + chain, + (pargs["umbrella-sampling"]) ? wf : WeightlessFunction() + ); + numeric_type = if pargs["numeric-type"] == "float64" + Float64; + elseif pargs["numeric-type"] == "float128" + Float128; + elseif pargs["numeric-type"] == "dec128" + Dec128; + elseif pargs["numeric-type"] == "big" + BigFloat; + else + error("numeric-type '$(pargs["numeric-type"])' not understood"); + exit(-1); + end + + # setup averagers + Avg, avgcons, avgconss = if pargs["umbrella-sampling"] + (UmbrellaAverager, + ((acc, chain) -> begin; + UmbrellaAverager( + StandardAverager( + numeric_type(0.0*acc(chain)), + numeric_type(0.0), + acc + ), + wf + ); + end), + ((acc, chain) -> begin; + UmbrellaAverager( + StandardAverager( + Vector{numeric_type}(0.0*acc(chain)), + numeric_type(0.0), + acc + ), + wf + ); + end) + ); + else + (StandardAverager, + ((acc, chain) -> begin; + StandardAverager( + numeric_type, + acc, chain + ) + end), + ((acc, chain) -> begin; + StandardAverager( + numeric_type, + Vector{numeric_type}, + acc, chain + ); + end) + ); + end + scalar_averagers = (Avg{numeric_type,numeric_type}[ + (avgcons(chain -> dot(end_to_end(chain), + end_to_end(chain)), chain)), + (avgcons(chain -> dot(chain_μ(chain), chain_μ(chain)), + chain)), + (avgcons(chain -> chain.U, chain)), + (avgcons(chain -> chain.U*chain.U, chain)), + ]); + vector_averagers = (Avg{Vector{numeric_type},numeric_type}[ + avgconss(end_to_end, chain), + avgconss(chain -> map(x -> x*x, end_to_end(chain)), chain), + avgconss(chain_μ, chain), + avgconss(chain -> map(x -> x*x, chain_μ(chain)), chain) + ]); + + outfile = open("$(pargs["prefix"])_trajectory.csv", "w"); + writedlm(outfile, ["step" "r1" "r2" "r3" "p1" "p2" "p3" "U"], ','); + rollfile = open("$(pargs["prefix"])_rolling.csv", "w"); + writedlm(rollfile, ["step" "r1" "r2" "r3" "r1sq" "r2sq" "r3sq" "rsq" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psq" "U" "Usq"], ','); + + start = time(); + last_update = start; + nacc = 0; + nacc_total = 0; + natt = 0; + + for step=1:nsteps + idx = rand(1:n(chain)); + dϕ = rand(dϕ_dist); + dθ = rand(dθ_dist); + trial_chain = EAPChain(chain); + successful = move!(trial_chain, idx, dϕ, dθ); + α = cluster_flip!(trial_chain, idx; ϵflip = pargs["cluster-prob"]); + if successful && acceptor(trial_chain, rand(); α = α) + chain = trial_chain; + nacc += 1; + nacc_total += 1; + end + natt += 1; + + if time() - last_update > pargs["update-freq"] + @info "elapsed: $(time() - start)"; + @info "step: $step / $nsteps"; + last_update = time(); + end + + if ( + pargs["step-adjust-scale"] != 1.0 && + step % pargs["steps-per-adjust"] == 0 + ) # adjust step size? + ar = nacc / natt; + if (ar > pargs["step-adjust-ub"] && + ϕstep != π && θstep != π/2) + @info "acceptance ratio is high; increasing step size"; + nacc = 0; + natt = 0; + ϕstep = min(π, ϕstep*pargs["step-adjust-scale"]); + θstep = min(π/2, θstep*pargs["step-adjust-scale"]); + elseif ar < pargs["step-adjust-lb"] + @info "acceptance ratio is low; decreasing step size"; + nacc = 0; + natt = 0; + ϕstep /= pargs["step-adjust-scale"]; + θstep /= pargs["step-adjust-scale"]; + end + dϕ_dist = Uniform(-ϕstep, ϕstep); + dθ_dist = Uniform(-θstep, θstep); + end + + foreach(a -> record!(a, chain), scalar_averagers); + foreach(a -> record!(a, chain), vector_averagers); + if step % pargs["stepout"] == 0 + writedlm(outfile, + hcat(step, transpose(chain.r), + transpose(chain_μ(chain)), chain.U), + ','); + writedlm(rollfile, + hcat( + step, + transpose(get_avg(vector_averagers[1])), # r + transpose(get_avg(vector_averagers[2])), # rj2 + get_avg(scalar_averagers[1]), # r2 + transpose(get_avg(vector_averagers[3])), # p + transpose(get_avg(vector_averagers[4])), # pj2 + get_avg(scalar_averagers[2]), # p2 + get_avg(scalar_averagers[3]), # U + get_avg(scalar_averagers[4]) # U2 + ), + ','); + + end + + end # steps + + AR = nacc_total / pargs["num-steps"]; + @info "total time elapsed: $(time() - start)"; + @info "acceptance rate: $AR"; + + close(outfile); + close(rollfile); + + return (scalar_averagers, vector_averagers, AR); +end + +(sas, vas, ar) = if pargs["profile"] + @info "Profiling the mcmc code..."; + mcmc(5, pargs); # run first to compile code + @profview mcmc(pargs["num-steps"], pargs); + println("Press ENTER to continue..."); + readline(stdin); + exit(1); +else + mcmc(pargs["num-steps"], pargs); +end + +println(" = $(get_avg(vas[1]))"); +println(" = $(get_avg(vas[1]) / (pargs["mlen"]*pargs["num-monomers"]))"); +println(" = $(get_avg(vas[2]))"); +println(" = $(get_avg(sas[1]))"); +println("

= $(get_avg(vas[3]))"); +println(" = $(get_avg(vas[4]))"); +println(" = $(get_avg(sas[2]))"); +println(" = $(get_avg(sas[3]))"); +println(" = $(get_avg(sas[4]))"); +println("AR = $ar"); diff --git a/mcmc_eap_chain.jl b/mcmc_eap_chain.jl index c6152e6..233abe4 100644 --- a/mcmc_eap_chain.jl +++ b/mcmc_eap_chain.jl @@ -11,6 +11,11 @@ include(joinpath("inc", "eap_chain.jl")); include(joinpath("inc", "average.jl")); include(joinpath("inc", "acceptance.jl")); +function flip_n!(ψs::Vector) + ψs[1] += π; + ψs[2] = π - ψs[2]; +end + s = ArgParseSettings(); @add_arg_table! s begin "--E0", "-e" @@ -85,7 +90,10 @@ s = ArgParseSettings(); arg_type = Float64; default = 3*π / 8; "--do-flips" - help = "try reflecting monomer about xy plane in trial moves" + help = "trial moves with flipping monomers" + action = :store_true + "--do-clustering" + help = "trial moves with clustering" action = :store_true "--theta-step", "-q" help = "maximum θ step length" From 4bdcab19710742ddfebdcfd70dc7130806e93dca Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Fri, 24 Sep 2021 13:57:08 -0400 Subject: [PATCH 02/17] updated README --- README.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/README.md b/README.md index 94bb30f..cc42a32 100644 --- a/README.md +++ b/README.md @@ -22,6 +22,12 @@ You can see options by running: julia mcmc_eap_chain.jl --help +An experimental algorithm for clustering of the chain is implmeneted in ``mcmc_clustering_eap_chain.jl``. +You can see options by running: + + julia mcmc_clustering_eap_chain.jl --help + + ## Questions Documentation of this work is still in progress. However, I am very receptive to questions and comments via email. You can email questions to grasingerm@gmail.com. From 24a7c07e9d5cb8e49e621955b4c80b48680f3512 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Fri, 24 Sep 2021 22:42:53 -0400 Subject: [PATCH 03/17] added some parameter studies; debugged adaptive step size --- mcmc_eap_chain.jl | 19 +++--- ...ting-compare-with-clustering_2021-09-24.jl | 62 +++++++++++++++++++ ...ninteracting-with-clustering_2021-09-24.jl | 42 +++++++++++++ 3 files changed, 116 insertions(+), 7 deletions(-) create mode 100644 run/noninteracting-compare-with-clustering_2021-09-24.jl create mode 100644 run/noninteracting-with-clustering_2021-09-24.jl diff --git a/mcmc_eap_chain.jl b/mcmc_eap_chain.jl index 233abe4..d4d0dd1 100644 --- a/mcmc_eap_chain.jl +++ b/mcmc_eap_chain.jl @@ -92,9 +92,6 @@ s = ArgParseSettings(); "--do-flips" help = "trial moves with flipping monomers" action = :store_true - "--do-clustering" - help = "trial moves with clustering" - action = :store_true "--theta-step", "-q" help = "maximum θ step length" arg_type = Float64; @@ -263,7 +260,9 @@ function mcmc(nsteps::Int, pargs) start = time(); last_update = start; - num_accepted = 0; + nacc = 0; + nacc_total = 0; + natt = 0; for init=1:pargs["num-inits"] if !force_ensemble_flag @@ -287,8 +286,10 @@ function mcmc(nsteps::Int, pargs) end if successful && acceptor(trial_chain, rand()) chain = trial_chain; - num_accepted += 1; + nacc += 1; + nacc_total += 1; end + natt += 1; if time() - last_update > pargs["update-freq"] @info "elapsed: $(time() - start)"; @@ -301,14 +302,18 @@ function mcmc(nsteps::Int, pargs) pargs["step-adjust-scale"] != 1.0 && step % pargs["steps-per-adjust"] == 0 ) # adjust step size? - acc_ratio = num_accepted / (step + (init - 1)*nsteps); + acc_ratio = nacc / natt; if (acc_ratio > pargs["step-adjust-ub"] && ϕstep != π && θstep != π/2) @info "acceptance ratio is high; increasing step size"; + nacc = 0; + natt = 0; ϕstep = min(π, ϕstep*pargs["step-adjust-scale"]); θstep = min(π/2, θstep*pargs["step-adjust-scale"]); elseif acc_ratio < pargs["step-adjust-lb"] @info "acceptance ratio is low; decreasing step size"; + nacc = 0; + natt = 0; ϕstep /= pargs["step-adjust-scale"]; θstep /= pargs["step-adjust-scale"]; end @@ -357,7 +362,7 @@ function mcmc(nsteps::Int, pargs) end - ar = num_accepted / (pargs["num-inits"]*pargs["num-steps"]); + ar = nacc_total / (pargs["num-inits"]*pargs["num-steps"]); @info "total time elapsed: $(time() - start)"; @info "acceptance rate: $ar"; diff --git a/run/noninteracting-compare-with-clustering_2021-09-24.jl b/run/noninteracting-compare-with-clustering_2021-09-24.jl new file mode 100644 index 0000000..3df16c1 --- /dev/null +++ b/run/noninteracting-compare-with-clustering_2021-09-24.jl @@ -0,0 +1,62 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = 0.0:0.5:10.0; +kTs = [1.0]; +Fzs = vcat(0.0:0.1:1.0, 1.5:0.5:5.0); +for Fz in Fzs, kT in kTs, E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => 1.0, :K2 => 0.0, + :kT => kT, :Fz => Fz, + :Fx => 0.0, :n => 100, :b => 1)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + + outfile = joinpath(workdir, "$(prefix(case))_standard.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case), "standard"))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "$(prefix(case))_umbrella.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting --umbrella-sampling -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)), "umbrella")`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "$(prefix(case))_clustering.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)), "clustering")`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); diff --git a/run/noninteracting-with-clustering_2021-09-24.jl b/run/noninteracting-with-clustering_2021-09-24.jl new file mode 100644 index 0000000..ee356d5 --- /dev/null +++ b/run/noninteracting-with-clustering_2021-09-24.jl @@ -0,0 +1,42 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = 0.0:1.0:5.0; +kTs = [1.0]; +Fzs = vcat(0.0:0.1:1.0, 1.5:0.5:5.0); +for Fz in Fzs, kT in kTs, E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => 1.0, :K2 => 0.0, + :kT => kT, :Fz => Fz, + :Fx => 0.0, :n => 100, :b => 1)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); From 2ab0aeff841a4017e556015daf773fa3194e0e8c Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Sat, 25 Sep 2021 22:04:00 -0400 Subject: [PATCH 04/17] debugged some parameter studies --- ...ting-compare-with-clustering_2021-09-24.jl | 72 +++++++++++++++++++ ...ting-compare-with-clustering_2021-09-24.jl | 11 ++- ...ninteracting-with-clustering_2021-09-24.jl | 2 + 3 files changed, 82 insertions(+), 3 deletions(-) create mode 100644 run/interacting-compare-with-clustering_2021-09-24.jl diff --git a/run/interacting-compare-with-clustering_2021-09-24.jl b/run/interacting-compare-with-clustering_2021-09-24.jl new file mode 100644 index 0000000..323308c --- /dev/null +++ b/run/interacting-compare-with-clustering_2021-09-24.jl @@ -0,0 +1,72 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))_run-$(fmt_int(case[:run]))"; +end + +cases = Any[]; +E0s = [1e-1; 1.0; 10.0]; +K1s = [1e-2; 1e-1; 1.0]; +kTs = [1.0; 0.1]; +Fzs = [0.0; 0.5; 1.0]; +Fxs = [0.0; 0.5; 1.0]; +ns = Int[100]; +bs = [0.5; 1.0; 2.0]; +runs = 1:25; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, K1 in K1s, run in runs + push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => 0.0, + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b, :run => run)); +end + +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); +mkpath(joinpath(workdir, "standard")); +mkpath(joinpath(workdir, "umbrella")); +mkpath(joinpath(workdir, "clustering")); + +pmap(case -> begin; + + outfile = joinpath(workdir, "$(prefix(case))_standard.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "standard", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "$(prefix(case))_umbrella.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type Ising --umbrella-sampling -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "umbrella", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "$(prefix(case))_clustering.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "clustering", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); diff --git a/run/noninteracting-compare-with-clustering_2021-09-24.jl b/run/noninteracting-compare-with-clustering_2021-09-24.jl index 3df16c1..7f48f10 100644 --- a/run/noninteracting-compare-with-clustering_2021-09-24.jl +++ b/run/noninteracting-compare-with-clustering_2021-09-24.jl @@ -27,9 +27,14 @@ end @info "total number of cases to run: $(length(cases))"; +mkpath(workdir); +mkpath(joinpath(workdir, "standard")); +mkpath(joinpath(workdir, "umbrella")); +mkpath(joinpath(workdir, "clustering")); + pmap(case -> begin; - outfile = joinpath(workdir, "$(prefix(case))_standard.out"); + outfile = joinpath(workdir, "standard", "$(prefix(case))_standard.out"); if !isfile(outfile) println("Running case: $case."); command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case), "standard"))`; @@ -39,7 +44,7 @@ pmap(case -> begin; println("Case: $case has already been run."); end - outfile = joinpath(workdir, "$(prefix(case))_umbrella.out"); + outfile = joinpath(workdir, "umbrella", "$(prefix(case))_umbrella.out"); if !isfile(outfile) println("Running case: $case."); command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting --umbrella-sampling -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)), "umbrella")`; @@ -49,7 +54,7 @@ pmap(case -> begin; println("Case: $case has already been run."); end - outfile = joinpath(workdir, "$(prefix(case))_clustering.out"); + outfile = joinpath(workdir, "clustering", "$(prefix(case))_clustering.out"); if !isfile(outfile) println("Running case: $case."); command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)), "clustering")`; diff --git a/run/noninteracting-with-clustering_2021-09-24.jl b/run/noninteracting-with-clustering_2021-09-24.jl index ee356d5..d9ab5ab 100644 --- a/run/noninteracting-with-clustering_2021-09-24.jl +++ b/run/noninteracting-with-clustering_2021-09-24.jl @@ -27,6 +27,8 @@ end @info "total number of cases to run: $(length(cases))"; +mkpath(workdir); + pmap(case -> begin; outfile = joinpath(workdir, "$(prefix(case)).out"); From 627a2b6745a3d3036ed111a3fb2145944fc84590 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Sun, 26 Sep 2021 21:09:15 -0400 Subject: [PATCH 05/17] tweaked two param studies --- ...ting-compare-with-clustering_2021-09-24.jl | 20 ++++--------------- ...ninteracting-with-clustering_2021-09-24.jl | 4 ++-- 2 files changed, 6 insertions(+), 18 deletions(-) diff --git a/run/noninteracting-compare-with-clustering_2021-09-24.jl b/run/noninteracting-compare-with-clustering_2021-09-24.jl index 7f48f10..5b50f39 100644 --- a/run/noninteracting-compare-with-clustering_2021-09-24.jl +++ b/run/noninteracting-compare-with-clustering_2021-09-24.jl @@ -16,28 +16,26 @@ end end cases = Any[]; -E0s = 0.0:0.5:10.0; +E0s = 0.0:1.0:5.0; kTs = [1.0]; -Fzs = vcat(0.0:0.1:1.0, 1.5:0.5:5.0); +Fzs = vcat(0.0:0.05:1.0, 1.5:0.5:5.0); for Fz in Fzs, kT in kTs, E0 in E0s push!(cases, Dict(:E0 => E0, :K1 => 1.0, :K2 => 0.0, :kT => kT, :Fz => Fz, :Fx => 0.0, :n => 100, :b => 1)); end - @info "total number of cases to run: $(length(cases))"; mkpath(workdir); mkpath(joinpath(workdir, "standard")); mkpath(joinpath(workdir, "umbrella")); -mkpath(joinpath(workdir, "clustering")); pmap(case -> begin; outfile = joinpath(workdir, "standard", "$(prefix(case))_standard.out"); if !isfile(outfile) println("Running case: $case."); - command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case), "standard"))`; + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "standard", prefix(case)))`; output = read(command, String); write(outfile, output); else @@ -47,17 +45,7 @@ pmap(case -> begin; outfile = joinpath(workdir, "umbrella", "$(prefix(case))_umbrella.out"); if !isfile(outfile) println("Running case: $case."); - command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting --umbrella-sampling -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)), "umbrella")`; - output = read(command, String); - write(outfile, output); - else - println("Case: $case has already been run."); - end - - outfile = joinpath(workdir, "clustering", "$(prefix(case))_clustering.out"); - if !isfile(outfile) - println("Running case: $case."); - command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)), "clustering")`; + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting --umbrella-sampling -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "umbrella", prefix(case)))`; output = read(command, String); write(outfile, output); else diff --git a/run/noninteracting-with-clustering_2021-09-24.jl b/run/noninteracting-with-clustering_2021-09-24.jl index d9ab5ab..12dd3a7 100644 --- a/run/noninteracting-with-clustering_2021-09-24.jl +++ b/run/noninteracting-with-clustering_2021-09-24.jl @@ -18,7 +18,7 @@ end cases = Any[]; E0s = 0.0:1.0:5.0; kTs = [1.0]; -Fzs = vcat(0.0:0.1:1.0, 1.5:0.5:5.0); +Fzs = vcat(0.0:0.05:1.0, 1.5:0.5:5.0); for Fz in Fzs, kT in kTs, E0 in E0s push!(cases, Dict(:E0 => E0, :K1 => 1.0, :K2 => 0.0, :kT => kT, :Fz => Fz, @@ -34,7 +34,7 @@ pmap(case -> begin; outfile = joinpath(workdir, "$(prefix(case)).out"); if !isfile(outfile) println("Running case: $case."); - command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 10000000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; output = read(command, String); write(outfile, output); else From d9fd77aaf91505a42d16c0ea681e03d6364f0820 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Mon, 27 Sep 2021 10:17:50 -0400 Subject: [PATCH 06/17] Update README.md --- README.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index cc42a32..a8bf114 100644 --- a/README.md +++ b/README.md @@ -26,7 +26,9 @@ An experimental algorithm for clustering of the chain is implmeneted in ``mcmc_c You can see options by running: julia mcmc_clustering_eap_chain.jl --help - + +It has been verified against analytical results, so I am relatively confident in its correctness. +The expectation is that the clustering algorithm will vastly improve the convergence rate for dielectric chains with monomer-monomer interactions, but this has not been studied yet. ## Questions Documentation of this work is still in progress. However, I am very receptive to questions and comments via email. You can email questions to grasingerm@gmail.com. From e6f93ca73fdac1e6ee3e18a3b215f259b03f564b Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Thu, 30 Sep 2021 13:09:55 -0400 Subject: [PATCH 07/17] added some cases to interacting-clustering study --- ...ting-compare-with-clustering_2021-09-28.jl | 108 ++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 run/interacting-compare-with-clustering_2021-09-28.jl diff --git a/run/interacting-compare-with-clustering_2021-09-28.jl b/run/interacting-compare-with-clustering_2021-09-28.jl new file mode 100644 index 0000000..ceb386e --- /dev/null +++ b/run/interacting-compare-with-clustering_2021-09-28.jl @@ -0,0 +1,108 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))_run-$(fmt_int(case[:run]))"; +end + +cases = Any[]; +#E0s = [1.0; 5.0]; +#K1s = [1e-1; 1.0]; +#kTs = [1.0; 0.1]; +#Fzs = [0.0; 0.5; 1.0]; +#Fxs = [0.0; 0.5; 1.0]; +#ns = Int[100]; +#bs = [1.0]; +runs = 1:25; +for run in runs + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 5, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 5, :K1 => 0.1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 0.1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 1, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 2, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 1, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 0.1, :Fz => 0, + :Fx => 1, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 0.1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 0.1, :Fz => 1, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 0.1, :K2 => 0.0, + :kT => 0.1, :Fz => 1, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 200, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 0.5, :run => run)); +end + +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); +mkpath(joinpath(workdir, "standard")); +mkpath(joinpath(workdir, "umbrella")); +mkpath(joinpath(workdir, "clustering")); + +pmap(case -> begin; + + outfile = joinpath(workdir, "standard", "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "standard", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "umbrella", "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type Ising --umbrella-sampling -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "umbrella", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "clustering", "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "clustering", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); From c936dfcf66efb1f70c50dd0558bc0de9be95b202 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Fri, 1 Oct 2021 15:08:45 -0400 Subject: [PATCH 08/17] added some cases to study --- ...ting-compare-with-clustering_2021-09-28.jl | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/run/interacting-compare-with-clustering_2021-09-28.jl b/run/interacting-compare-with-clustering_2021-09-28.jl index ceb386e..5810932 100644 --- a/run/interacting-compare-with-clustering_2021-09-28.jl +++ b/run/interacting-compare-with-clustering_2021-09-28.jl @@ -64,6 +64,36 @@ for run in runs push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, :kT => 1, :Fz => 0, :Fx => 0, :n => 100, :b => 0.5, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 2.0, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 1, + :Fx => 0, :n => 100, :b => 0.5, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 1, :n => 100, :b => 0.5, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 1, + :Fx => 1, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 3, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 3, + :Fx => 3, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 3, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 5, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 3, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 5, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 5, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 5, + :Fx => 5, :n => 100, :b => 1, :run => run)); end @info "total number of cases to run: $(length(cases))"; From 5e9dedd02ff0017d8231f35b5e260ee5e345869f Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Tue, 26 Oct 2021 11:36:38 -0400 Subject: [PATCH 09/17] added an ad hoc and experimental 2D version of the mcmc code for verification with analytical results --- 2D/inc/acceptance.jl | 39 ++ 2D/inc/algorithm_help_msg.jl | 28 ++ 2D/inc/average.jl | 124 +++++ 2D/inc/dielectric_dipole.jl | 4 + 2D/inc/dipole_response.jl | 27 ++ 2D/inc/eap_chain.jl | 278 +++++++++++ 2D/inc/energy.jl | 23 + 2D/inc/free_energy.jl | 25 + 2D/inc/interacting_energy.jl | 4 + 2D/inc/langevin.jl | 23 + 2D/inc/ni_energy.jl | 3 + 2D/inc/output.jl | 22 + 2D/inc/solvers.jl | 430 ++++++++++++++++++ 2D/inc/types.jl | 4 + 2D/mcmc_clustering_eap_chain.jl | 322 +++++++++++++ ...ting-compare-with-clustering_2021-09-24.jl | 72 +++ ...ting-compare-with-clustering_2021-09-28.jl | 138 ++++++ ...interacting-dielectric-study_2021-10-24.jl | 115 +++++ 2D/run/interacting_dielectric_study.jl | 46 ++ 2D/run/interacting_polar_study.jl | 45 ++ ...ting-compare-with-clustering_2021-09-24.jl | 55 +++ ...ninteracting-with-clustering_2021-09-24.jl | 44 ++ 2D/run/prelim/callbacks.jl | 31 ++ 2D/run/prelim/param_study.jl | 44 ++ 2D/run/prelim/param_study2.jl | 47 ++ 2D/run/prelim/param_study3.jl | 44 ++ 2D/run/prelim/param_study4.jl | 44 ++ 2D/run/prelim/param_study5.jl | 45 ++ 2D/run/prelim/param_study6.jl | 45 ++ 2D/run/prelim/param_study7.jl | 45 ++ 2D/run/prelim/param_study8.jl | 45 ++ 2D/run/prelim/param_study9.jl | 44 ++ 2D/run/prelim/test_noE.jl | 44 ++ mcmc_clustering_eap_chain.jl | 2 +- 34 files changed, 2350 insertions(+), 1 deletion(-) create mode 100644 2D/inc/acceptance.jl create mode 100644 2D/inc/algorithm_help_msg.jl create mode 100644 2D/inc/average.jl create mode 100644 2D/inc/dielectric_dipole.jl create mode 100644 2D/inc/dipole_response.jl create mode 100644 2D/inc/eap_chain.jl create mode 100644 2D/inc/energy.jl create mode 100644 2D/inc/free_energy.jl create mode 100644 2D/inc/interacting_energy.jl create mode 100644 2D/inc/langevin.jl create mode 100644 2D/inc/ni_energy.jl create mode 100644 2D/inc/output.jl create mode 100644 2D/inc/solvers.jl create mode 100644 2D/inc/types.jl create mode 100644 2D/mcmc_clustering_eap_chain.jl create mode 100644 2D/run/interacting-compare-with-clustering_2021-09-24.jl create mode 100644 2D/run/interacting-compare-with-clustering_2021-09-28.jl create mode 100644 2D/run/interacting-dielectric-study_2021-10-24.jl create mode 100644 2D/run/interacting_dielectric_study.jl create mode 100644 2D/run/interacting_polar_study.jl create mode 100644 2D/run/noninteracting-compare-with-clustering_2021-09-24.jl create mode 100644 2D/run/noninteracting-with-clustering_2021-09-24.jl create mode 100644 2D/run/prelim/callbacks.jl create mode 100644 2D/run/prelim/param_study.jl create mode 100644 2D/run/prelim/param_study2.jl create mode 100644 2D/run/prelim/param_study3.jl create mode 100644 2D/run/prelim/param_study4.jl create mode 100644 2D/run/prelim/param_study5.jl create mode 100644 2D/run/prelim/param_study6.jl create mode 100644 2D/run/prelim/param_study7.jl create mode 100644 2D/run/prelim/param_study8.jl create mode 100644 2D/run/prelim/param_study9.jl create mode 100644 2D/run/prelim/test_noE.jl diff --git a/2D/inc/acceptance.jl b/2D/inc/acceptance.jl new file mode 100644 index 0000000..63285b4 --- /dev/null +++ b/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 diff --git a/2D/inc/algorithm_help_msg.jl b/2D/inc/algorithm_help_msg.jl new file mode 100644 index 0000000..688e7c4 --- /dev/null +++ b/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 :) +"""; diff --git a/2D/inc/average.jl b/2D/inc/average.jl new file mode 100644 index 0000000..4334153 --- /dev/null +++ b/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 diff --git a/2D/inc/dielectric_dipole.jl b/2D/inc/dielectric_dipole.jl new file mode 100644 index 0000000..d78581b --- /dev/null +++ b/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 diff --git a/2D/inc/dipole_response.jl b/2D/inc/dipole_response.jl new file mode 100644 index 0000000..88b7765 --- /dev/null +++ b/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 diff --git a/2D/inc/eap_chain.jl b/2D/inc/eap_chain.jl new file mode 100644 index 0000000..e193275 --- /dev/null +++ b/2D/inc/eap_chain.jl @@ -0,0 +1,278 @@ +using Distributions; +using NLsolve; +using ForwardDiff; + +include(joinpath(@__DIR__, "dipole_response.jl")); + +ϕ_dist = Uniform(0.0, 2*π); + +# the way that the energy functions are organized is a mess; I'm sorry +# also... not sorry +abstract type Energy end + +mutable struct EAPChain + b::Float64; + E0::Float64; + μ::DipoleResponse; + UFunction::Energy; + kT::Float64; + Fz::Float64; + Fx::Float64; + ϕs::FdVector; + cϕs::FdVector; + sϕs::FdVector; + n̂s::FdMatrix; + μs::FdMatrix; + us::FdVector; + xs::FdMatrix; + r::FdVector; + U::Float64; +end + +@inline n(chain::EAPChain) = length(chain.ϕs); + +@inline n̂(cϕ, sϕ) = [cϕ; sϕ]; + +@inline n̂j(chain::EAPChain, idx::Int) = n̂(chain.cϕs[idx], chain.sϕs[idx]); + +function update_xs!(chain::EAPChain, idx::Int) + @warn "This implementation of update_xs! is generally slower than the vectorized implementation"; + if idx == 1 + @inbounds chain.xs[:, idx] = chain.b / 2 * n̂j(chain, idx); + idx += 1; + end + for i=idx:n(chain) + @inbounds chain.xs[:, i] = (chain.xs[:, i-1] + + chain.b / 2 * (n̂j(chain, i) + n̂j(chain, i-1))); + end +end + +@inline function update_xs!(chain::EAPChain) + chain.xs[:, :] = chain.b*(cumsum(chain.n̂s, dims=2) - 0.5*chain.n̂s); +end + +#= +function update_xs!(chain::EAPChain) + @show chain.n̂s + @show cumsum(chain.n̂s, dims=2); + @show (cumsum(chain.n̂s, dims=2) - 0.5*chain.n̂s); + idx = rand(1:n(chain)); + @assert(dot(chain.n̂s[:, idx], chain.n̂s[:, idx]) ≈ 1.0); + chain.xs[:, :] = chain.b*(cumsum(chain.n̂s, dims=2) - 0.5*chain.n̂s); +end +=# + +@inline u(E0::Real, μ::FdVector) = -1/2*E0*μ[1]; + +function EAPChain(pargs::Dict) + ϕs = rand(ϕ_dist, pargs["num-monomers"]); + + dr = if pargs["chain-type"] == "dielectric" + DielectricResponse(pargs["K1"], pargs["K2"]); + elseif pargs["chain-type"] == "polar" + PolarResponse(pargs["mu"] * [1.0 0.0; 0.0 1.0]); + else + error("chain-type is not understood."); + end + + ret = EAPChain( + pargs["mlen"], + pargs["E0"], + dr, + if pargs["energy-type"] == "noninteracting" + NonInteractingEnergy(); + elseif pargs["energy-type"] == "interacting" + InteractingEnergy(); + elseif pargs["energy-type"] == "Ising" + IsingEnergy(); + else + error("energy-type is not understood."); + end, + pargs["kT"], + pargs["Fz"], + pargs["Fx"], + ϕs, + map(cos, ϕs), + map(sin, ϕs), + zeros(2, pargs["num-monomers"]), + zeros(2, pargs["num-monomers"]), + zeros(pargs["num-monomers"]), + zeros(2, pargs["num-monomers"]), + zeros(2), + 0.0 + ); + ret.n̂s[:, :] = hcat(map(j -> n̂j(ret, j), 1:n(ret))...); + @simd for i=1:n(ret) + @inbounds ret.μs[:, i] = dr(ret.E0, ret.cϕs[i], ret.sϕs[i]); + end + ret.us[:] = map(i -> u(ret.E0, view(ret.μs, :, i)), 1:n(ret)); + update_xs!(ret); + ret.r[:] = end_to_end(ret); + ret.U = U(ret); + return ret; +end + +function EAPChain(chain::EAPChain) + return EAPChain( + chain.b, + chain.E0, + chain.μ, + chain.UFunction, + chain.kT, + chain.Fz, + chain.Fx, + chain.ϕs[:], + chain.cϕs[:], + chain.sϕs[:], + chain.n̂s[:, :], + chain.μs[:, :], + chain.us[:], + chain.xs[:, :], + chain.r[:], + chain.U + ); +end + +# TODO: consider rewriting this to make better use of past calculations +# this is probably by far the slowest calculation *shrug emoji* +function U_interaction(chain::EAPChain) + U = 0.0; + @inbounds for i=1:n(chain) + @simd for j=i+1:n(chain) # once debugged, use inbounds + r = chain.xs[:, i] - chain.xs[:, j]; + r2 = dot(r, r); + rmag = sqrt(r2); + r̂ = r / rmag; + r3 = r2*rmag; + μi = view(chain.μs, :, i); + μj = view(chain.μs, :, j); + U += (dot(μi, μj) - 3*dot(μi, r̂)*dot(μj, r̂)) / (4*π*r3); + end + end + return U; +end + +# TODO: consider rewriting this to make better use of past calculations +# this is probably by far the slowest calculation *shrug emoji* +function U_Ising(chain::EAPChain) + U = 0.0; + @inbounds for i=1:n(chain)-1 + r = chain.xs[:, i] - chain.xs[:, i+1]; + r2 = dot(r, r); + rmag = sqrt(r2); + r̂ = r / rmag; + r3 = r2*rmag; + μi = view(chain.μs, :, i); + μj = view(chain.μs, :, i+1); + U += (dot(μi, μj) - 3*dot(μi, r̂)*dot(μj, r̂)) / (4*π*r3); + end + return U; +end + +function move!(chain::EAPChain, idx::Int, dϕ::Real) + @inbounds begin + chain.ϕs[idx] += dϕ; + chain.cϕs[idx] = cos(chain.ϕs[idx]); + chain.sϕs[idx] = sin(chain.ϕs[idx]); + + # the order of these updates is important + chain.n̂s[:, idx] = n̂j(chain, idx); + chain.μs[:, idx] = chain.μ(chain.E0, chain.cϕs[idx], chain.sϕs[idx]); + chain.us[idx] = u(chain.E0, view(chain.μs, :, idx)); + update_xs!(chain); + chain.r[:] = end_to_end(chain); + chain.U = U(chain); + end + return true; +end + +function flip_n!(chain::EAPChain, idx::Int) + move!(chain, idx, π) +end + +pflip_linear(ip::Real) = (1 + ip) / 2; + +function cluster_flip!(chain::EAPChain, idx::Int; + pflip::Function = pflip_linear, + ϵflip::Real = 1.0, + ηerr::Real = 1e-10 + ) + + # grow the cluster to the right + upper_p = NaN; + upper_idx = idx; + while true + if upper_idx >= n(chain) + upper_p = 0.0; + break; + end + upper_p = pflip_linear(dot(n̂j(chain, upper_idx), n̂j(chain, upper_idx+1))); + if rand() <= upper_p + upper_idx += 1 + else + break; + end + end + @assert(!isnan(upper_p) && upper_idx <= n(chain)); + + # grow the cluster to the left + lower_p = NaN; + lower_idx = idx; + while true + if lower_idx <= 1 + lower_p = 0.0; + break; + end + lower_p = pflip_linear(dot(n̂j(chain, lower_idx), n̂j(chain, lower_idx-1))); + if rand() <= lower_p + lower_idx -= 1 + else + break; + end + end + @assert(!isnan(lower_p) && lower_idx >= 1); + + if rand() <= ϵflip + prev_n̂s = copy(chain.n̂s[:, lower_idx:upper_idx]); + for i in lower_idx:upper_idx + flip_n!(chain, i); + end + + if norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) > ηerr + println("prev"); + display(prev_n̂s) + println(); + println("curr"); + display(chain.n̂s[:, lower_idx:upper_idx]) + println(); + @show(norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf)) + end + @assert(norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) < ηerr); + + new_upper_p = if upper_idx < n(chain) + pflip_linear(dot(n̂j(chain, upper_idx), n̂j(chain, upper_idx+1))); + else + 0 + end + new_lower_p = if lower_idx > 1 + pflip_linear(dot(n̂j(chain, lower_idx), n̂j(chain, lower_idx-1))); + else + 0 + end + return ( + ((1 - new_upper_p)*(1 - new_lower_p)) / + ((1 - upper_p)*(1 - lower_p)) + ); + else + return 1.0; + end + +end + +@inline end_to_end(chain::EAPChain) = @inbounds (chain.xs[:, end] + + chain.b/2.0*n̂j(chain, n(chain))); + +@inline chain_μ(chain::EAPChain) = reshape(sum(chain.μs, dims=2), 2); + +include(joinpath(@__DIR__, "energy.jl")); +@inline U(chain::EAPChain) = chain.UFunction(chain); # forward energy call to chain diff --git a/2D/inc/energy.jl b/2D/inc/energy.jl new file mode 100644 index 0000000..b87798d --- /dev/null +++ b/2D/inc/energy.jl @@ -0,0 +1,23 @@ +include(joinpath(@__DIR__, "types.jl")); + +struct NonInteractingEnergy <: Energy end + +(::Energy)(::EAPChain) = error("Not yet implemented."); + +@inline function (::NonInteractingEnergy)(chain::EAPChain) + return (sum(chain.us) - dot(end_to_end(chain), [chain.Fx; chain.Fz])); +end + +struct InteractingEnergy <: Energy end + +@inline function (::InteractingEnergy)(chain::EAPChain) + return (sum(chain.us) + U_interaction(chain) + - dot(end_to_end(chain), [chain.Fx; chain.Fz])); +end + +struct IsingEnergy <: Energy end + +@inline function (::IsingEnergy)(chain::EAPChain) + return (sum(chain.us) + U_Ising(chain) + - dot(end_to_end(chain), [chain.Fx; chain.Fz])); +end diff --git a/2D/inc/free_energy.jl b/2D/inc/free_energy.jl new file mode 100644 index 0000000..68e85f9 --- /dev/null +++ b/2D/inc/free_energy.jl @@ -0,0 +1,25 @@ +function free_energy(xs::Vector) + I1 = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*u(θs), + _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + I2 = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*log(ρ(θs, xs)), + _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + return I1[1] + kT*I2[1] - kT*N*log(N); +end + +function free_energy_smallomega(x::Vector) + c, λ, α = x[:]; + I1 = pcubature(θs -> c*(1 - ω(θs) + α*cos(θs[1])*sin(θs[2]))*exp(λ*cos(θs[2]))*sin(θs[2])*u(θs), + _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + I2 = pcubature(θs -> c*(1 - ω(θs) + α*cos(θs[1])*sin(θs[2]))*exp(λ*cos(θs[2]))*sin(θs[2])*log(ρ(θs, xs)), + _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + return I1[1] + kT*I2[1] - kT*N*log(N); +end + +function free_energy_smalllambda(x::Vector) + c, λ, α = x[:]; + I1 = pcubature(θs -> c*(1 + λ*cos(θs[2]) + α*cos(θs[1])*sin(θs[2]))*exp(-ω(θs))*sin(θs[2])*u(θs), + _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + I2 = pcubature(θs -> c*(1 + λ*cos(θs[2]) + α*cos(θs[1])*sin(θs[2]))*exp(-ω(θs))*sin(θs[2])*log(ρ(θs, xs)), + _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + return I1[1] + kT*I2[1] - kT*N*log(N); +end diff --git a/2D/inc/interacting_energy.jl b/2D/inc/interacting_energy.jl new file mode 100644 index 0000000..72fad07 --- /dev/null +++ b/2D/inc/interacting_energy.jl @@ -0,0 +1,4 @@ +@inline function U(chain::EAPChain) + return (sum(chain.us) + U_interaction(chain) + - dot(end_to_end(chain), [chain.Fx; chain.Fz])); +end diff --git a/2D/inc/langevin.jl b/2D/inc/langevin.jl new file mode 100644 index 0000000..50a8efa --- /dev/null +++ b/2D/inc/langevin.jl @@ -0,0 +1,23 @@ +function langevin_inv(y::Real) + return (3*y - y/5 * (6*y^2 + y^4 - 2*y^6))/(1 - y^2); +end + +function free_energy_langevin(k::Real, T::Real, N::Int, γ::Real) + linvγ = langevin_inv(γ); + return N*k*T*( log( linvγ / (4*pi*sinh(linvγ)) ) + γ*linvγ ); +end + +function force_langevin(k::Real, T::Real, b::Real, γ::Real) + return k*T/b * langevin_inv(γ); +end + +function forces(As::Vector, γs::Vector, N::Int) + n = length(As) + grads = zeros(n); + grads[1] = (As[2]-As[1])/(γs[2]-γs[1]); + grads[end] = (As[n]-As[n-1])/(γs[n]-γs[n-1]); + for i=2:(n-1) + grads[i] = (As[i+1] - As[i-1])/(γs[i+1] - γs[i-1]); + end + return grads / N; +end diff --git a/2D/inc/ni_energy.jl b/2D/inc/ni_energy.jl new file mode 100644 index 0000000..ce1efce --- /dev/null +++ b/2D/inc/ni_energy.jl @@ -0,0 +1,3 @@ +@inline function U(chain::EAPChain) + return (sum(chain.us) - dot(chain.r, [chain.Fx; 0.0; chain.Fz])); +end diff --git a/2D/inc/output.jl b/2D/inc/output.jl new file mode 100644 index 0000000..48b93f8 --- /dev/null +++ b/2D/inc/output.jl @@ -0,0 +1,22 @@ +include("free_energy.jl"); + +function standard_solution_output(xs::Vector) + println(); + println("A = $(free_energy(xs))"); + println("lambda = $(xs[2])"); + println("alpha = $(xs[3])"); + println("C = $(xs[1])"); + println("|lambda| = $(sqrt(xs[2]^2 + xs[3]^2))"); + μ = polarization(xs); + println("mu1 = $(μ[1])"); + println("mu2 = $(μ[2])"); + println("mu3 = $(μ[3])"); +end + +function residuals_output(rs::Real) + println("L2 = $rs"); +end + +function residuals_output(rs::Vector) + println("L2 = $(norm(rs, 2))"); +end diff --git a/2D/inc/solvers.jl b/2D/inc/solvers.jl new file mode 100644 index 0000000..4a30cd0 --- /dev/null +++ b/2D/inc/solvers.jl @@ -0,0 +1,430 @@ +include("langevin.jl"); + +using NLopt; +using Cubature; +using LinearAlgebra; + + _ALGO_DICT_ = Dict("PRAXIS" => :LN_PRAXIS, "SBPLX" => :LN_SBPLX, + "ISRES" => :GN_ISRES, "Newton" => :NEWTON, + "SLSQP" => :LD_SLSQP, "LBFGS" => :LD_LBFGS, + "AUGLAG" => :AUGLAG, "Anneal" => :ANNEAL); + +function algostr_to_symbol(algostr::AbstractString) + if algostr in keys(_ALGO_DICT_) + return _ALGO_DICT_[algostr]; + else + error("Do not understand the algorithm: \"" * algostr * "\""); + end + return :AJKLLKJALDSKJFAS +end + +function parse_algostr(algostr::AbstractString) + if startswith(algostr, "try") + steps = map(strip, split(algostr[4:end], "then")); + return map(parse_algostr, steps); + end + return map(algostr_to_symbol, map(strip, split(algostr, ","))); +end + +function parse_max_iters(max_iters_str::AbstractString) + if startswith(max_iters_str, "try") + steps = map(strip, split(max_iters_str[4:end], "then")); + return map(parse_max_iters, steps); + end + return map(x -> parse(Int, x), map(strip, split(max_iters_str, ","))); +end + +function initial_guess(γx::Real, γz::Real) + γ = sqrt(γx^2 + γz^2); + return [(γ > 0) ? N*langevin_inv(γ)*csch(langevin_inv(γ))/(4*pi) : 0; + (γz > 0) ? γz / γ * langevin_inv(γ) : 0; + (γx > 0) ? γx / γ * langevin_inv(γ) : 0]; +end + +function initial_guess(γ::Real) + return [γ > 0 ? N*langevin_inv(γ)*csch(langevin_inv(γ))/(4*pi) : 0; + γ > 0 ? langevin_inv(γ) : 0; + 0]; +end + +function initial_guess(ω::Real, γs::Tuple{Real,Real}) + sqrtω = sqrt(abs(ω)); + erfsω = ω > 0 ? erf(sqrtω) : erfi(sqrtω); + expω = exp(ω); + sqrtpi = sqrt(pi); + γx, γz = γs; + return [ + N * sqrtω / (2*pi*sqrtpi*erfsω); + (2*expω*sqrtpi*γz*ω*erfω) / (-2*sqrtω + expω*sqrtpi*erfsω); + (4*expω*sqrtpi*γx*ω*erfω) / (2*sqrtω + expω*sqrtpi*(-1+2*ω)*erfsω) + ]; +end + +function initial_guess(ωs::Tuple{Real,Real}, γ::Real) + ωx, ωz = ωs; + ω = ωx + ωz; + sqrtωx = sqrt(abs(ωx)); + sqrtωz = sqrt(abs(ωz)); + sqrtω = sqrt(abs(ω)); + erfsω = ω > 0 ? erf(sqrtω) : erfi(sqrtω); + expω = exp(ω); + sqrtpi = sqrt(pi); + return [ + N * sqrtω / (2*pi*sqrtpi*erfsω); + (2*expω*sqrtpi*sqrtωz*sqrtω*erfsω) / (-2*sqrtω + expω*sqrtpi*erfsω); + (4*expω*sqrtpi*sqrtωx*sqrtω*erfsω) / (2*sqrtω + expω*sqrtpi*(-1+2*ω)*erfsω) + ]; +end + + _RHS_ECOORD_ = Float64[N; γz*N; γx*N]; + +function residuals(Is::Vector) + global _RHS_; + map(i -> Is[i][1] - _RHS_[i], 1:3) +end + +ω(θs::Vector) = ω(θs[1], θs[2]); + +ρ(ϕ::Real, θ::Real, c::Real, λ::Real, α::Real) = c * exp(-ω(ϕ, θ) + λ*cos(θ) + + α*cos(ϕ)*sin(θ)); +ρ(θs::Vector, xs::Vector) = ρ(θs[1], θs[2], xs[1], xs[2], xs[3]); +u(ϕ::Real, θ::Real) = kT * ω(ϕ, θ) + _u0_; +u(θs::Vector) = kT * ω(θs[1], θs[2]) + _u0_; + +integrand_functions(xs::Vector) = ( + [ + (θs::Vector) -> ρ(θs, xs)*sin(θs[2]), + (θs::Vector) -> ρ(θs, xs)*sin(θs[2])*cos(θs[2]), + (θs::Vector) -> ρ(θs, xs)*sin(θs[2])*sin(θs[2])*cos(θs[1]) + ] +); + + _LOWER_BOUNDS_ = [0.0, 0.0]; + _UPPER_BOUNDS_ = [2*pi, pi]; + +function integrate_sys(xs::Vector) + return map(integrand -> pcubature(integrand, _LOWER_BOUNDS_, _UPPER_BOUNDS_; + reltol=RELTOL_INT, maxevals=MAXEVALS_INT), + integrand_functions(xs)); +end + +function gradc(xs::Vector) + global RELTOL_INT; + global MAXEVALS_INT; + + grad = zeros(3); + grad[1], = pcubature(θs -> ρ(θs, xs)/xs[1]*sin(θs[2]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + grad[2], = pcubature(θs -> ρ(θs, xs)/xs[1]*sin(θs[2])*cos(θs[2]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + grad[3], = pcubature(θs -> ρ(θs, xs)/xs[1]*sin(θs[2])*sin(θs[2])*cos(θs[1]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + + return grad; +end + +function gradλ(xs::Vector) + global RELTOL_INT; + global MAXEVALS_INT; + + grad = zeros(3); + grad[1], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*cos(θs[2]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + grad[2], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*cos(θs[2])*cos(θs[2]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + grad[3], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*sin(θs[2])*cos(θs[1])*cos(θs[2]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + + return grad; +end + +function gradα(xs::Vector) + global RELTOL_INT; + global MAXEVALS_INT; + + grad = zeros(3); + grad[1], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*sin(θs[2])*cos(θs[1]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + grad[2], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*cos(θs[2])*sin(θs[2])*cos(θs[1]), _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + grad[3], = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*(sin(θs[2])*cos(θs[1]))^2, _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT); + + return grad; +end + +function hessian(xs::Vector) + h = zeros(3, 3); + + h[:, 1] = gradc(xs); + h[:, 2] = gradλ(xs); + h[:, 3] = gradα(xs); + + return h; +end + +function newton_raphson_iteration(xs::Vector, rs::Vector) + h = hessian(xs); + return xs - inv(h) * rs; +end + +function newton_raphson(xs::Vector; reltol::Real=1e-10, abstol::Real=1e-8, max_iters::Int=100000) + initial_Is = integrate_sys(xs); + current_residuals = residuals(initial_Is); + residual_norm = norm(current_residuals); + init_residual_norm = residual_norm; + + if residual_norm < abstol + return (xs, current_residuals, true, 0); + end + + global print_every_percent_compl; + iters_per_msg = ceil(Int, print_every_percent_compl * max_iters); + + for iter = 1:max_iters + + lin_alg_error_flag = false; + try + xs = newton_raphson_iteration(xs, current_residuals); + catch + lin_alg_error_flag = true; + end + + # if Hessian is singular, or Newton iteration fails, perturb guess + if lin_alg_error_flag + @warn "hessian was singular; perturbing guess" + xs += rand(3); + end + + Is = integrate_sys(xs); + current_residuals = residuals(Is); + residual_norm = norm(current_residuals); + + relnorm = residual_norm / init_residual_norm; + + if residual_norm < abstol || relnorm < reltol + return (xs, current_residuals, true, iter); + end + + if iter % iters_per_msg == 0 + println("iter: $iter / $max_iters; Rnorm = $residual_norm"); + end + + end + + return (xs, current_residuals, false, max_iters); +end + +function simulated_annealing(xs::Vector; anneal_factor::Real=1e-3, + init_T::Real=10.0, final_T::Real=1e-7, + relax_steps::Int=100, + init_delta::Real=0.5, + delta_adapt_factor::Real=1.2, + adapt_acc_thrshld::Real=0.7, + reltol::Real=1e-10, abstol::Real=1e-8, + max_iters::Int=1000000) + + initial_Is = integrate_sys(xs); + current_residuals = residuals(initial_Is); + residual_norm = norm(current_residuals); + init_residual_norm = residual_norm; + Is = initial_Is; + + if residual_norm < abstol + return (xs, current_residuals, true, 0); + end + + acc = 0; + delta = init_delta; + T = init_T; + + for iter = 1:max_iters + + # check to see if relaxation should occur + if iter % relax_steps == 0 + T *= (1.0 - anneal_factor); + if T < final_T + break; + end + + # adapt step size to current acceptance ratio + if acc / relax_steps > adapt_acc_thrshld + delta *= delta_adapt_factor; + elseif acc / relax_steps < (1 - adapt_acc_thrshld) + delta /= delta_adapt_factor; + end + + # reset acceptance count + acc = 0; + end + + # trial move + choice = rand(1:3); + dx = (rand(Float64)*2 - 1) * delta; + new_xs = copy(xs); + new_xs[choice] += dx; + + new_Is = integrate_sys(new_xs); + new_residuals = residuals(new_Is); + + # metropolis acceptance criteria + if rand(Float64) < exp(-(dot(new_residuals, new_residuals) - + dot(current_residuals, current_residuals)) / T) + + # accept move + acc += 1; + copy!(xs, new_xs); + current_residuals = new_residuals; + Is = copy(new_Is); + + # check for convergence + residual_norm = norm(current_residuals); + relnorm = residual_norm / init_residual_norm; + if residual_norm < abstol || relnorm < reltol + return (xs, current_residuals, true, iter); + end + + end + + end + + return (xs, current_residuals, false, iter); + +end + +_GRADS_ = [gradc, gradλ, gradα]; + +NLopt_raint_functions = map(i -> begin + (xs::Vector, grad::Vector) -> begin + if length(grad) > 0 + grad = gradc(xs); + end + + return (pcubature(integrand_functions(xs)[i], _LOWER_BOUNDS_, _UPPER_BOUNDS_; + reltol=RELTOL_INT, maxevals=MAXEVALS_INT)[1] - _RHS_[i]); + end +end, 1:3); + +function NLopt_optimization_function!(xs::Vector, grad::Vector) + + Is = integrate_sys(xs); + current_residuals = residuals(Is); + + if length(grad) > 0 + h = hessian(xs); + grad = vec(transpose(current_residuals) * h); + end + + return 0.5 * dot(current_residuals, current_residuals); +end + + _CONV_SMBLS_ = [:SUCCESS; :STOPVAL_REACHED; :FTOL_REACHED; :XTOL_REACHED]; + +function run_solver(algo::Symbol; + x0s::Vector = DEFAULT_INIT_GUESS, reltol::Real=1e-10, + abstol::Real=1e-8, max_iters::Int=10000) + solver = init_solver(algo, x0s, reltol, abstol, max_iters); + return solve(solver); +end + +function run_solver(algos::Vector{Symbol}, maxs_iters::Vector{Int}; + x0s::Vector = DEFAULT_INIT_GUESS, reltol::Real=1e-10, + abstol::Real=1e-8) + + if length(algos) != length(maxs_iters) + error("Length of algorithm and maximum iteration vectors must be equal"); + end + + total_iterations = 0; + xs, rs, converged = x0s, [NaN; NaN; NaN], false; + + for (algo, max_iters) in zip(algos, maxs_iters) + (xs, rs, converged, iterations) = run_solver(algo; x0s=xs, + reltol=reltol, + abstol=abstol, + max_iters=max_iters); + total_iterations += iterations; + end + + return (xs, rs, converged, total_iterations); +end + +function run_solver(algo_steps::Vector{Vector{Symbol}}, + max_iters_by_step::Vector{Vector{Int}}; + x0s::Vector = DEFAULT_INIT_GUESS, reltol::Real=1e-10, + abstol::Real=1e-6) + + if length(algo_steps) != length(max_iters_by_step) + error("Length of algorithm and maximum iteration vectors must be equal"); + end + + total_iterations = 0; + xs, rs, converged = x0s, [NaN; NaN; NaN], false; + + for (algo_step, max_iters_step) in zip(algo_steps, max_iters_by_step) + (xs, rs, converged, iterations) = run_solver(algo_step, max_iters_step; + x0s=x0s, + reltol=reltol, + abstol=abstol); + total_iterations += iterations; + if converged + return (xs, rs, converged, total_iterations); + end + end + + return (xs, rs, converged, total_iterations); +end + + _NLOPT_SOLVERS_ = [:LN_PRAXIS; :LN_SBPLX; :GN_ISRES; :LD_SLSQP; + :LD_LBFGS; :AUGLAG]; + _NLOPT_REQ_BOUNDS_ = [:GN_ISRES]; + _NLOPT_REQ_RAINTS_ = [:GN_ISRES; :LD_SLSQP; :AUGLAG]; + _NLOPT_DEFAULT_LOWER_BOUNDS_ = [-N; -100.0; -100.0]; + _NLOPT_DEFAULT_UPPER_BOUNDS_ = [100.0*N; 100.0; 100.0]; + +function init_solver(algo::Symbol, x0s::Vector, reltol::Real, abstol::Real, + max_iters::Int) + + global _NUMEVALS_; + global _ITERS_PER_PRINT_; + global print_every_percent_compl; + + _NUMEVALS_ = 0; + _ITERS_PER_PRINT_ = ceil(Int, print_every_percent_compl * max_iters); + + if algo in _NLOPT_SOLVERS_ + opt = Opt(algo, 3); + min_objective!(opt, NLopt_optimization_function!); + if algo in _NLOPT_REQ_BOUNDS_ + lower_bounds!(_NLOPT_DEFAULT_LOWER_BOUNDS_); + upper_bounds!(_NLOPT_DEFAULT_UPPER_BOUNDS_); + end + if algo in _NLOPT_REQ_RAINTS_ + for i=1:3 + equality_raint!(opt, NLopt_raint_functions[i]); + end + end + if algo == :AUGLAG + local_opt = init_solver(:LD_LBFGS, x0s, reltol, abstol, max_iters); + local_optimizer!(opt, local_opt[1]); + end + ftol_rel!(opt, reltol); + ftol_abs!(opt, abstol); + maxeval!(opt, max_iters); + return (opt, x0s); + elseif algo == :NEWTON + return ((xs) -> newton_raphson(xs; reltol=reltol, abstol=abstol, max_iters=max_iters), + x0s); + elseif algo == :ANNEAL + return ((xs) -> simulated_annealing(xs; anneal_factor=ANNEAL_FACTOR, + init_T=ANNEAL_INIT_T, final_T=ANNEAL_FINAL_T, + relax_steps=ANNEAL_RELAX_STEPS, + init_delta=ANNEAL_INIT_DELTA, + delta_adapt_factor=ANNEAL_ADAPT_FACTOR, + adapt_acc_thrshld=ANNEAL_ACC_THRSHLD, + reltol=reltol, abstol=abstol, + max_iters=max_iters), + x0s); + else + error("Do not understand algorithm symbol: ", algo); + end +end + +function solve(opt_pair::Tuple{NLopt.Opt, Vector}) + (minf, minx, ret) = optimize(opt_pair[1], opt_pair[2]); + if ret == :XTOL_REACHED; @warn "Convergence given by x-tolerance"; end + return (minx, sqrt(2 * minf), ret in _CONV_SMBLS_, _NUMEVALS_); +end + +solve(newton_pair::Tuple{Function, Vector}) = newton_pair[1](newton_pair[2]); diff --git a/2D/inc/types.jl b/2D/inc/types.jl new file mode 100644 index 0000000..5e8df39 --- /dev/null +++ b/2D/inc/types.jl @@ -0,0 +1,4 @@ +CnMatrix = AbstractMatrix{Int}; +FdMatrix = AbstractMatrix{Float64}; +CnVector = AbstractVector{Int}; +FdVector = AbstractVector{Float64}; diff --git a/2D/mcmc_clustering_eap_chain.jl b/2D/mcmc_clustering_eap_chain.jl new file mode 100644 index 0000000..b3f0e7d --- /dev/null +++ b/2D/mcmc_clustering_eap_chain.jl @@ -0,0 +1,322 @@ +using ArgParse; +using Distributions; +using LinearAlgebra; +using Logging; +using DelimitedFiles; +#using ProfileView; +using DecFP; +using Quadmath; + +include(joinpath("inc", "eap_chain.jl")); +include(joinpath("inc", "average.jl")); +include(joinpath("inc", "acceptance.jl")); + +s = ArgParseSettings(); +@add_arg_table! s begin + "--E0", "-e" + help = "magnitude of electric field" + arg_type = Float64 + default = 0.0 + "--chain-type", "-T" + help = "chain type (dielectric|polar)" + arg_type = String + default = "dielectric" + "--K1", "-J" + help = "dipole susceptibility along the monomer axis (dielectric chain)" + arg_type = Float64 + default = 1.0 + "--K2", "-K" + help = "dipole susceptibility orthogonal to the monomer axis (dielectric chain)" + arg_type = Float64 + default = 0.0 + "--mu", "-m" + arg_type = Float64 + default = 1e-2 + help = "dipole magnitude (electret chain)" + "--energy-type", "-u" + help = "energy type (noninteracting|interacting|Ising)" + arg_type = String + default = "noninteracting" + "--kT", "-k" + help = "dimensionless temperature" + arg_type = Float64 + default = 1.0 + "--Fz", "-F" + help = "force in the z-direction (direction of E-field; force ensemble)" + arg_type = Float64 + default = 0.0 + "--Fx", "-G" + help = "force in the x-direction (force ensemble)" + arg_type = Float64 + default = 0.0 + "--mlen", "-b" + help = "monomer length" + arg_type = Float64 + default = 1.0 + "--num-monomers", "-n" + help = "number of monomers" + arg_type = Int + default = 100 + "--num-steps", "-N" + help = "number of steps" + arg_type = Int + default = convert(Int, 1e6) + "--phi-step", "-p" + help = "maximum ϕ step length" + arg_type = Float64; + default = 3*π / 8; +# "--do-clustering" +# help = "trial moves with clustering" +# action = :store_true + "--cluster-prob" + help = "probability of flipping a cluster" + arg_type = Float64; + default = 0.5; + "--step-adjust-lb", "-L" + help = "adjust step sizes if acc. ratio below this threshold" + arg_type = Float64 + default = 0.15 + "--step-adjust-ub", "-U" + help = "adjust step sizes if acc. ratio above this threshold" + arg_type = Float64 + default = 0.40 + "--step-adjust-scale", "-A" + help = "scale factor for adjusting step sizes (> 1.0)" + arg_type = Float64 + default = 1.1 + "--steps-per-adjust", "-S" + help = "steps between step size adjustments" + arg_type = Int + default = 2500 + "--umbrella-sampling", "-B" + help = "use umbrella sampling (w/ electrostatic weight function)" + action = :store_true + "--update-freq" + help = "update frequency (seconds)" + arg_type = Float64; + default = 15.0; + "--verbose", "-v" + help = "verbosity level: 0-nothing, 1-errors, 2-warnings, 3-info" + arg_type = Int + default = 3 + "--prefix", "-P" + help = "prefix for output files" + arg_type = String + default = "eap-mcmc" + "--postfix", "-Q" + help = "postfix for output files" + arg_type = String + default = "" + "--stepout", "-s" + help = "steps between storing microstates" + arg_type = Int + default = 500 + "--numeric-type" + help = "numerical data type for averaging (float64|float128|dec128|big)" + arg_type = String + default = "float64" + "--profile", "-Z" + help = "profile the program" + action = :store_true +end + +pargs = parse_args(s); + +if pargs["verbose"] == 3 + global_logger(ConsoleLogger(stderr, Logging.Info)); +elseif pargs["verbose"] == 2 + global_logger(ConsoleLogger(stderr, Logging.Warn)); +elseif pargs["verbose"] == 1 + global_logger(ConsoleLogger(stderr, Logging.Error)); +else + global_logger(Logging.NullLogger()); +end + +function mcmc(nsteps::Int, pargs) + ϕstep = pargs["phi-step"]; + dϕ_dist = Uniform(-ϕstep, ϕstep); + chain = EAPChain(pargs); + chain.U = U(chain); + wf = AntiDipoleWeightFunction(chain); + acceptor = Metropolis( + chain, + (pargs["umbrella-sampling"]) ? wf : WeightlessFunction() + ); + numeric_type = if pargs["numeric-type"] == "float64" + Float64; + elseif pargs["numeric-type"] == "float128" + Float128; + elseif pargs["numeric-type"] == "dec128" + Dec128; + elseif pargs["numeric-type"] == "big" + BigFloat; + else + error("numeric-type '$(pargs["numeric-type"])' not understood"); + exit(-1); + end + + # setup averagers + Avg, avgcons, avgconss = if pargs["umbrella-sampling"] + (UmbrellaAverager, + ((acc, chain) -> begin; + UmbrellaAverager( + StandardAverager( + numeric_type(0.0*acc(chain)), + numeric_type(0.0), + acc + ), + wf + ); + end), + ((acc, chain) -> begin; + UmbrellaAverager( + StandardAverager( + Vector{numeric_type}(0.0*acc(chain)), + numeric_type(0.0), + acc + ), + wf + ); + end) + ); + else + (StandardAverager, + ((acc, chain) -> begin; + StandardAverager( + numeric_type, + acc, chain + ) + end), + ((acc, chain) -> begin; + StandardAverager( + numeric_type, + Vector{numeric_type}, + acc, chain + ); + end) + ); + end + scalar_averagers = (Avg{numeric_type,numeric_type}[ + (avgcons(chain -> dot(end_to_end(chain), + end_to_end(chain)), chain)), + (avgcons(chain -> dot(chain_μ(chain), chain_μ(chain)), + chain)), + (avgcons(chain -> chain.U, chain)), + (avgcons(chain -> chain.U*chain.U, chain)), + ]); + vector_averagers = (Avg{Vector{numeric_type},numeric_type}[ + avgconss(end_to_end, chain), + avgconss(chain -> map(x -> x*x, end_to_end(chain)), chain), + avgconss(chain_μ, chain), + avgconss(chain -> map(x -> x*x, chain_μ(chain)), chain) + ]); + + outfile = open("$(pargs["prefix"])_trajectory.csv", "w"); + writedlm(outfile, ["step" "r1" "r3" "p1" "p3" "U"], ','); + rollfile = open("$(pargs["prefix"])_rolling.csv", "w"); + writedlm(rollfile, ["step" "r1" "r3" "r1sq" "r3sq" "rsq" "p1" "p3" "p1sq" "p3sq" "psq" "U" "Usq"], ','); + + start = time(); + last_update = start; + nacc = 0; + nacc_total = 0; + natt = 0; + + for step=1:nsteps + idx = rand(1:n(chain)); + dϕ = rand(dϕ_dist); + trial_chain = EAPChain(chain); + successful = move!(trial_chain, idx, dϕ); + α = cluster_flip!(trial_chain, idx; ϵflip = pargs["cluster-prob"]); + if successful && acceptor(trial_chain, rand(); α = α) + chain = trial_chain; + nacc += 1; + nacc_total += 1; + end + natt += 1; + + if time() - last_update > pargs["update-freq"] + @info "elapsed: $(time() - start)"; + @info "step: $step / $nsteps"; + last_update = time(); + end + + if ( + pargs["step-adjust-scale"] != 1.0 && + step % pargs["steps-per-adjust"] == 0 + ) # adjust step size? + ar = nacc / natt; + if (ar > pargs["step-adjust-ub"] && + ϕstep != π) + @info "acceptance ratio is high; increasing step size"; + nacc = 0; + natt = 0; + ϕstep = min(π, ϕstep*pargs["step-adjust-scale"]); + elseif ar < pargs["step-adjust-lb"] + @info "acceptance ratio is low; decreasing step size"; + nacc = 0; + natt = 0; + ϕstep /= pargs["step-adjust-scale"]; + end + dϕ_dist = Uniform(-ϕstep, ϕstep); + end + + foreach(a -> record!(a, chain), scalar_averagers); + foreach(a -> record!(a, chain), vector_averagers); + if step % pargs["stepout"] == 0 + writedlm(outfile, + hcat(step, transpose(chain.r), + transpose(chain_μ(chain)), chain.U), + ','); + writedlm(rollfile, + hcat( + step, + transpose(get_avg(vector_averagers[1])), # r + transpose(get_avg(vector_averagers[2])), # rj2 + get_avg(scalar_averagers[1]), # r2 + transpose(get_avg(vector_averagers[3])), # p + transpose(get_avg(vector_averagers[4])), # pj2 + get_avg(scalar_averagers[2]), # p2 + get_avg(scalar_averagers[3]), # U + get_avg(scalar_averagers[4]) # U2 + ), + ','); + + end + + end # steps + + AR = nacc_total / pargs["num-steps"]; + @info "total time elapsed: $(time() - start)"; + @info "acceptance rate: $AR"; + + close(outfile); + close(rollfile); + + return (scalar_averagers, vector_averagers, AR); +end + +#= +(sas, vas, ar) = if pargs["profile"] + @info "Profiling the mcmc code..."; + mcmc(5, pargs); # run first to compile code + @profview mcmc(pargs["num-steps"], pargs); + println("Press ENTER to continue..."); + readline(stdin); + exit(1); +else + mcmc(pargs["num-steps"], pargs); +end +=# +(sas, vas, ar) = mcmc(pargs["num-steps"], pargs); + +println(" = $(get_avg(vas[1]))"); +println(" = $(get_avg(vas[1]) / (pargs["mlen"]*pargs["num-monomers"]))"); +println(" = $(get_avg(vas[2]))"); +println(" = $(get_avg(sas[1]))"); +println("

= $(get_avg(vas[3]))"); +println(" = $(get_avg(vas[4]))"); +println(" = $(get_avg(sas[2]))"); +println(" = $(get_avg(sas[3]))"); +println(" = $(get_avg(sas[4]))"); +println("AR = $ar"); diff --git a/2D/run/interacting-compare-with-clustering_2021-09-24.jl b/2D/run/interacting-compare-with-clustering_2021-09-24.jl new file mode 100644 index 0000000..323308c --- /dev/null +++ b/2D/run/interacting-compare-with-clustering_2021-09-24.jl @@ -0,0 +1,72 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))_run-$(fmt_int(case[:run]))"; +end + +cases = Any[]; +E0s = [1e-1; 1.0; 10.0]; +K1s = [1e-2; 1e-1; 1.0]; +kTs = [1.0; 0.1]; +Fzs = [0.0; 0.5; 1.0]; +Fxs = [0.0; 0.5; 1.0]; +ns = Int[100]; +bs = [0.5; 1.0; 2.0]; +runs = 1:25; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, K1 in K1s, run in runs + push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => 0.0, + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b, :run => run)); +end + +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); +mkpath(joinpath(workdir, "standard")); +mkpath(joinpath(workdir, "umbrella")); +mkpath(joinpath(workdir, "clustering")); + +pmap(case -> begin; + + outfile = joinpath(workdir, "$(prefix(case))_standard.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "standard", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "$(prefix(case))_umbrella.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type Ising --umbrella-sampling -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "umbrella", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "$(prefix(case))_clustering.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "clustering", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); diff --git a/2D/run/interacting-compare-with-clustering_2021-09-28.jl b/2D/run/interacting-compare-with-clustering_2021-09-28.jl new file mode 100644 index 0000000..5810932 --- /dev/null +++ b/2D/run/interacting-compare-with-clustering_2021-09-28.jl @@ -0,0 +1,138 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))_run-$(fmt_int(case[:run]))"; +end + +cases = Any[]; +#E0s = [1.0; 5.0]; +#K1s = [1e-1; 1.0]; +#kTs = [1.0; 0.1]; +#Fzs = [0.0; 0.5; 1.0]; +#Fxs = [0.0; 0.5; 1.0]; +#ns = Int[100]; +#bs = [1.0]; +runs = 1:25; +for run in runs + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 5, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 5, :K1 => 0.1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 0.1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 1, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 2, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 1, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 0.1, :Fz => 0, + :Fx => 1, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 0.1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 0.1, :Fz => 1, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 0.1, :K2 => 0.0, + :kT => 0.1, :Fz => 1, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 200, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 0.5, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 2.0, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 1, + :Fx => 0, :n => 100, :b => 0.5, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 1, :n => 100, :b => 0.5, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 1, + :Fx => 1, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 3, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 3, + :Fx => 3, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 3, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 5, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 3, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 5, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 5, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 5, + :Fx => 5, :n => 100, :b => 1, :run => run)); +end + +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); +mkpath(joinpath(workdir, "standard")); +mkpath(joinpath(workdir, "umbrella")); +mkpath(joinpath(workdir, "clustering")); + +pmap(case -> begin; + + outfile = joinpath(workdir, "standard", "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "standard", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "umbrella", "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type Ising --umbrella-sampling -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "umbrella", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "clustering", "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "clustering", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); diff --git a/2D/run/interacting-dielectric-study_2021-10-24.jl b/2D/run/interacting-dielectric-study_2021-10-24.jl new file mode 100644 index 0000000..47af575 --- /dev/null +++ b/2D/run/interacting-dielectric-study_2021-10-24.jl @@ -0,0 +1,115 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))_run-$(fmt_int(case[:run]))"; +end + +cases = Any[]; +#E0s = [1.0; 5.0]; +#K1s = [1e-1; 1.0]; +#kTs = [1.0; 0.1]; +#Fzs = [0.0; 0.5; 1.0]; +#Fxs = [0.0; 0.5; 1.0]; +#ns = Int[100]; +#bs = [1.0]; +runs = 1:25; +for run in runs + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 5, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 5, :K1 => 0.1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 0.1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 1, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 2, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 1, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 0.1, :Fz => 0, + :Fx => 1, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 0.1, :Fz => 0, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 0.1, :Fz => 1, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 0.1, :K2 => 0.0, + :kT => 0.1, :Fz => 1, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 200, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 0.5, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 0, :n => 100, :b => 2.0, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 1, + :Fx => 0, :n => 100, :b => 0.5, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 1, :n => 100, :b => 0.5, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 1, + :Fx => 1, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 3, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 3, + :Fx => 3, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 3, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 5, + :Fx => 0, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 3, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 1, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 0, + :Fx => 5, :n => 100, :b => 1, :run => run)); + push!(cases, Dict(:E0 => 5, :K1 => 1, :K2 => 0.0, + :kT => 1, :Fz => 5, + :Fx => 5, :n => 100, :b => 1, :run => run)); +end + +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); + +pmap(case -> begin; + + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); diff --git a/2D/run/interacting_dielectric_study.jl b/2D/run/interacting_dielectric_study.jl new file mode 100644 index 0000000..cd534cd --- /dev/null +++ b/2D/run/interacting_dielectric_study.jl @@ -0,0 +1,46 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [1e-1; 1.0; 10.0]; +K1s = [0.0; 1e-1; 0.5; 1.0; 2.0]; +K2s = [0.0; 1e-1; 0.5; 1.0; 2.0]; +kTs = [1.0]; +Fzs = [0.0; 0.5; 1.0; 2.0]; +Fxs = [0.0; 0.5; 1.0; 2.0]; +ns = Int[100; 200]; +bs = [0.5; 1.0; 2.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, K1 in K1s, K2 in K2s + if K1==K2; continue; end + push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2, + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type interacting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 250000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end +end, cases); diff --git a/2D/run/interacting_polar_study.jl b/2D/run/interacting_polar_study.jl new file mode 100644 index 0000000..009bd22 --- /dev/null +++ b/2D/run/interacting_polar_study.jl @@ -0,0 +1,45 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_mu-$(fmt(case[:mu]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [0.0; 1e-1; 1.0; 10.0]; +μs = [1e-6; 1e-3; 1e-2; 1e-1; 0.5; 1.0; 2.0; 10.0]; +kTs = [1.0]; +Fzs = sort(vcat(-[0.25; 0.5; 0.75; 1.0; 2.0; 3.0; 4.0; 5.0; 10.0; 20.0], + [0.0; 0.25; 0.5; 0.75; 1.0; 2.0; 3.0; 4.0; 5.0; 10.0; 20.0])); +Fxs = [0.0; 0.25; 0.5; 0.75; 1.0; 2.0; 3.0; 4.0; 5.0; 10.0; 20.0]; +ns = Int[100; 200]; +bs = [0.5; 1.0; 2.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, μ in μs + push!(cases, Dict(:E0 => E0, :mu => μ, + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type polar --energy-type interacting -b $(case[:b]) --E0 $(case[:E0]) --mu $(case[:mu]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case))) --numeric-type big`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end +end, cases); diff --git a/2D/run/noninteracting-compare-with-clustering_2021-09-24.jl b/2D/run/noninteracting-compare-with-clustering_2021-09-24.jl new file mode 100644 index 0000000..5b50f39 --- /dev/null +++ b/2D/run/noninteracting-compare-with-clustering_2021-09-24.jl @@ -0,0 +1,55 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = 0.0:1.0:5.0; +kTs = [1.0]; +Fzs = vcat(0.0:0.05:1.0, 1.5:0.5:5.0); +for Fz in Fzs, kT in kTs, E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => 1.0, :K2 => 0.0, + :kT => kT, :Fz => Fz, + :Fx => 0.0, :n => 100, :b => 1)); +end +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); +mkpath(joinpath(workdir, "standard")); +mkpath(joinpath(workdir, "umbrella")); + +pmap(case -> begin; + + outfile = joinpath(workdir, "standard", "$(prefix(case))_standard.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "standard", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + + outfile = joinpath(workdir, "umbrella", "$(prefix(case))_umbrella.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting --umbrella-sampling -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, "umbrella", prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); diff --git a/2D/run/noninteracting-with-clustering_2021-09-24.jl b/2D/run/noninteracting-with-clustering_2021-09-24.jl new file mode 100644 index 0000000..12dd3a7 --- /dev/null +++ b/2D/run/noninteracting-with-clustering_2021-09-24.jl @@ -0,0 +1,44 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = 0.0:1.0:5.0; +kTs = [1.0]; +Fzs = vcat(0.0:0.05:1.0, 1.5:0.5:5.0); +for Fz in Fzs, kT in kTs, E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => 1.0, :K2 => 0.0, + :kT => kT, :Fz => Fz, + :Fx => 0.0, :n => 100, :b => 1)); +end + +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); + +pmap(case -> begin; + + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 10000000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); diff --git a/2D/run/prelim/callbacks.jl b/2D/run/prelim/callbacks.jl new file mode 100644 index 0000000..fef4a4e --- /dev/null +++ b/2D/run/prelim/callbacks.jl @@ -0,0 +1,31 @@ +include(joinpath(@__DIR__, "..", "eap_chain.jl")); + +const stepout = 500; +outpath(pargs) = "$(pargs["prefix"])_trajectory.csv"; + +[ + (chain::EAPChain, step::Int, startup::Bool, cleanup::Bool, pargs) -> begin; + global outfile; + if startup + outfile = open(outpath(pargs), "a"); + end + end, + + (chain::EAPChain, step::Int, startup::Bool, cleanup::Bool, pargs) -> begin; + global outfile; + if cleanup + close(outfile); + end + end, + + (chain::EAPChain, step::Int, startup::Bool, cleanup::Bool, pargs) -> begin; + global outfile; + if step % stepout == 0 + writedlm(outfile, + hcat(step, transpose(chain.r), + transpose(chain_μ(chain)), chain.U), + ','); + end + end + +] diff --git a/2D/run/prelim/param_study.jl b/2D/run/prelim/param_study.jl new file mode 100644 index 0000000..6bbf3f6 --- /dev/null +++ b/2D/run/prelim/param_study.jl @@ -0,0 +1,44 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [10.0; 1.0; 0.1; 0.01; 1e-3; 0.0]; +K1s = [10.0; 1.0; 0.1; 0.0]; +K2s = [10.0; 1.0; 0.1; 0.0]; +kTs = [100.0; 10.0; 1.0; 0.1; 1e-3]; +Fzs = [0.0; 1e-3; 1.0; 10.0]; +Fxs = [0.0; 1e-3; 1.0; 10.0]; +ns = Int[10; 100; 250; 1000]; +bs = [1.0; 0.1; 10.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, K2 in K2s, K1 in K1s, + E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2, :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 interacting_chain.jl -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end +end, cases); diff --git a/2D/run/prelim/param_study2.jl b/2D/run/prelim/param_study2.jl new file mode 100644 index 0000000..8e050bc --- /dev/null +++ b/2D/run/prelim/param_study2.jl @@ -0,0 +1,47 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [10.0; 1.0; 1e-3]; +Ks = zip([1.0; 0.0; 0.5; 1.0], [0.0; 1.0; 1.0; 0.5]) +kTs = [100.0; 1.0; 1e-3]; +Fzs = [0.0; 1e-3; 1.0; 10.0]; +Fxs = [0.0; 1e-3; 1.0; 10.0]; +ns = Int[10; 100; 250; 1000]; +bs = [1.0; 0.1; 10.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, Kpair in Ks, + E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => Kpair[1], :K2 => Kpair[2], + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + for iter=1:10 + outfile = joinpath(workdir, "$(prefix(case))_iter-$iter.out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 interacting_chain.jl -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 100000 -v 2 --prefix $(joinpath(workdir, prefix(case)*"_iter-$iter"))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case, $iter has already been run."); + end + end +end, cases); diff --git a/2D/run/prelim/param_study3.jl b/2D/run/prelim/param_study3.jl new file mode 100644 index 0000000..c8dc93a --- /dev/null +++ b/2D/run/prelim/param_study3.jl @@ -0,0 +1,44 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [0.01; 0.05; 0.1; 0.5; 1.0; 2.0; 5.0; 10.0; 20.0]; +Ks = zip([1.0, 0.0], [0.0, 1.0]) +kTs = [1.0;]; +Fzs = [0.0; 0.025; 0.05; 0.075; 0.1; 0.5; 1.0; 2.5; 5.0; 10.0; 20.0]; +Fxs = [0.0; 0.025; 0.05; 0.075; 0.1; 0.5; 1.0; 2.5; 5.0; 10.0; 20.0]; +ns = Int[100]; +bs = [1.0;]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, Kpair in Ks, + E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => Kpair[1], :K2 => Kpair[2], + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 ni_chain.jl -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 500000 --num-inits 1 --step-adjust-scale 1.05 --step-adjust-ub 0.39 --step-adjust-lb 0.19 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end +end, cases); diff --git a/2D/run/prelim/param_study4.jl b/2D/run/prelim/param_study4.jl new file mode 100644 index 0000000..ae9d518 --- /dev/null +++ b/2D/run/prelim/param_study4.jl @@ -0,0 +1,44 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [0.0]; +Ks = zip([1.0], [0.0]) +kTs = [1.0;]; +Fzs = vcat(collect(range(0.0, 1.0; step=0.05)), range(2.0, 100.0; step=2.0)); +Fxs = [0.0;]; +ns = Int[100]; +bs = [1.0;]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, Kpair in Ks, + E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => Kpair[1], :K2 => Kpair[2], + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 ni_chain.jl -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 250000 --num-inits 1 --step-adjust-scale 1.05 --step-adjust-ub 0.39 --step-adjust-lb 0.19 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end +end, cases); diff --git a/2D/run/prelim/param_study5.jl b/2D/run/prelim/param_study5.jl new file mode 100644 index 0000000..6cff403 --- /dev/null +++ b/2D/run/prelim/param_study5.jl @@ -0,0 +1,45 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [1.0; 0.1; 1e-3]; +Ks = zip([1.0; 0.0], [0.0; 1.0]) +kTs = [1.0]; +Fzs = vcat(collect(range(0.0, 1.0; step=0.1)), range(5.0, 100.0; step=5.0)); +Fxs = [0.0]; +ns = Int[100]; +bs = [1.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, Kpair in Ks, + E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => Kpair[1], :K2 => Kpair[2], + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 250000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case, $iter has already been run."); + end +end, cases); diff --git a/2D/run/prelim/param_study6.jl b/2D/run/prelim/param_study6.jl new file mode 100644 index 0000000..26f9f0e --- /dev/null +++ b/2D/run/prelim/param_study6.jl @@ -0,0 +1,45 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = range(1.5, 2.5; length=100); +Ks = zip([1.0; 0.0], [0.0; 1.0]) +kTs = [1.0]; +Fzs = vcat(collect(range(0.0, 1.0; step=0.5)), range(5.0, 10.0; step=2.5)); +Fxs = [0.0]; +ns = Int[100]; +bs = [1.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, Kpair in Ks, + E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => Kpair[1], :K2 => Kpair[2], + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 250000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end +end, cases); diff --git a/2D/run/prelim/param_study7.jl b/2D/run/prelim/param_study7.jl new file mode 100644 index 0000000..6a0ff83 --- /dev/null +++ b/2D/run/prelim/param_study7.jl @@ -0,0 +1,45 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [1.0; 1.25; 1.5; 2.0; 2.5; 5.0; 7.5; 10.0]; +Ks = zip([1.0; 0.0], [0.0; 1.0]); +kTs = [1.0]; +Fzs = vcat(collect(range(0.0, 1.0; step=0.1)), range(2.0, 4.0; step=1.0), range(5.0, 20.0; step=2.5), 40.0); +Fxs = [0.0]; +ns = Int[100]; +bs = [1.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, Kpair in Ks, + E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => Kpair[1], :K2 => Kpair[2], + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting --do-flips -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1500000 -v 2 --prefix $(joinpath(workdir, prefix(case))) --umbrella-sampling --numeric-type big`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end +end, cases); diff --git a/2D/run/prelim/param_study8.jl b/2D/run/prelim/param_study8.jl new file mode 100644 index 0000000..9f8999b --- /dev/null +++ b/2D/run/prelim/param_study8.jl @@ -0,0 +1,45 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = vcat(range(sqrt(2)-0.75, sqrt(2)-0.25; length=10), range(sqrt(2)-0.25, sqrt(2)+0.25; length=25), range(sqrt(2)+0.25, sqrt(2)+0.75; length=10)); +Ks = zip([1.0; 0.0], [0.0; 1.0]); +kTs = [1.0]; +Fzs = [5.0; 2.5; 1.0]; +Fxs = [0.0]; +ns = Int[100]; +bs = [1.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, Kpair in Ks, + E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => Kpair[1], :K2 => Kpair[2], + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting --do-flips -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 2000000 -v 2 --prefix $(joinpath(workdir, prefix(case))) --umbrella-sampling --numeric-type big`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end +end, cases); diff --git a/2D/run/prelim/param_study9.jl b/2D/run/prelim/param_study9.jl new file mode 100644 index 0000000..3004179 --- /dev/null +++ b/2D/run/prelim/param_study9.jl @@ -0,0 +1,44 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = vcat(range(0.0, 0.6; length=10), range(sqrt(2)-0.75, sqrt(2)-0.25; length=10), range(sqrt(2)-0.25, sqrt(2)+0.25; length=25), range(sqrt(2)+0.25, sqrt(2)+0.75; length=10), range(2.5, 5.0; length=20)); +#Ks = zip([0.0; 1.0]); +kTs = [1.0]; +Fzs = [5.0; 10.0; 25.0; 50.0; 100.0]; +Fxs = [0.0]; +ns = Int[100]; +bs = [1.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => 0.0, :K2 => 1.0, + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_eap_chain.jl --chain-type dielectric --energy-type noninteracting -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 2500000 -v 2 --prefix $(joinpath(workdir, prefix(case))) --numeric-type big`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end +end, cases); diff --git a/2D/run/prelim/test_noE.jl b/2D/run/prelim/test_noE.jl new file mode 100644 index 0000000..1a383fa --- /dev/null +++ b/2D/run/prelim/test_noE.jl @@ -0,0 +1,44 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [0.0]; +Ks = zip([1.0], [0.0]) +kTs = [1.0;]; +Fzs = vcat(collect(range(0.0, 1.0; step=0.05)), range(2.0, 100.0; step=2.0)); +Fxs = [0.0;]; +ns = Int[100]; +bs = [1.0;]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, Kpair in Ks, + E0 in E0s + push!(cases, Dict(:E0 => E0, :K1 => Kpair[1], :K2 => Kpair[2], + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b)); +end + +@info "total number of cases to run: $(length(cases))"; + +pmap(case -> begin; + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 ni_chain.jl -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 10000000 --num-inits 1 --step-adjust-scale 1.1 --step-adjust-ub 0.39 --step-adjust-lb 0.19 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end +end, cases); diff --git a/mcmc_clustering_eap_chain.jl b/mcmc_clustering_eap_chain.jl index 430f0fb..85e0803 100644 --- a/mcmc_clustering_eap_chain.jl +++ b/mcmc_clustering_eap_chain.jl @@ -3,7 +3,7 @@ using Distributions; using LinearAlgebra; using Logging; using DelimitedFiles; -using ProfileView; +#using ProfileView; using DecFP; using Quadmath; From 96ff38b8900395bea8350076019d3cc2eb2e7682 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Wed, 27 Oct 2021 09:56:51 -0400 Subject: [PATCH 10/17] tweaked directory structure and README --- README.md | 54 ++++++++++--------- .../aggregate_mcmc.jl | 0 scripts/compare_mcmc.jl | 40 ++++++++++++++ .../create_graph_mean_field.jl | 0 .../install_dependencies.jl | 0 scripts/statistical_analysis_mcmc.jl | 40 ++++++++++++++ 6 files changed, 108 insertions(+), 26 deletions(-) rename aggregate_mcmc.jl => scripts/aggregate_mcmc.jl (100%) create mode 100644 scripts/compare_mcmc.jl rename create_graph_mean_field.jl => scripts/create_graph_mean_field.jl (100%) rename install_dependencies.jl => scripts/install_dependencies.jl (100%) create mode 100644 scripts/statistical_analysis_mcmc.jl diff --git a/README.md b/README.md index a8bf114..2924df3 100644 --- a/README.md +++ b/README.md @@ -6,9 +6,9 @@ As of now, chains consisting of polar monomers and dipole-dipole interactions (f ## Getting started This research code is written in the [julia language](https://julialang.org/). -In the root directory, there is a script for installing package dependences: +In the ``scripts`` directory, there is a script for installing package dependences: - julia install_dependencies.jl + julia scripts/install_dependencies.jl The main file for the mean-field theory simulations is ``mean_field.jl``. You can see options by running: @@ -35,36 +35,38 @@ Documentation of this work is still in progress. However, I am very receptive to ## Citations There are a few papers associated with polymer-stats. -If you found this code useful for your research, please be kind enough to cite it, using the DOIs 10.1039/D0SM00845A and 10.1073/pnas.2102477, or the following BiBTeX entries: +If you found this code useful for your research, please be kind enough to cite it, using the DOIs 10.1039/D0SM00845A, 10.1073/pnas.2102477 and 10.1016/j.jmps.2021.104658, or the following BiBTeX entries: @article{grasinger2020statistical, - title={Statistical mechanical analysis of the electromechanical coupling in an electrically-responsive polymer chain}, - author={Grasinger, Matthew and Dayal, Kaushik}, - journal={Soft Matter}, - volume={16}, - number={27}, - pages={6265--6284}, - year={2020}, - publisher={Royal Society of Chemistry}, - doi={10.1039/D0SM00845A} + title={Statistical mechanical analysis of the electromechanical coupling in an electrically-responsive polymer chain}, + author={Grasinger, Matthew and Dayal, Kaushik}, + journal={Soft Matter}, + volume={16}, + number={27}, + pages={6265--6284}, + year={2020}, + publisher={Royal Society of Chemistry}, + doi={10.1039/D0SM00845A} } @article {grasinger2021flexoelectricity, - author = {Grasinger, Matthew and Mozaffari, Kosar and Sharma, Pradeep}, - title = {Flexoelectricity in soft elastomers and the molecular mechanisms underpinning the design and emergence of giant flexoelectricity}, - volume = {118}, - number = {21}, - year = {2021}, - doi = {10.1073/pnas.2102477118}, - publisher = {National Academy of Sciences}, - issn = {0027-8424}, - journal = {Proceedings of the National Academy of Sciences} + author = {Grasinger, Matthew and Mozaffari, Kosar and Sharma, Pradeep}, + title = {Flexoelectricity in soft elastomers and the molecular mechanisms underpinning the design and emergence of giant flexoelectricity}, + volume = {118}, + number = {21}, + year = {2021}, + doi = {10.1073/pnas.2102477118}, + publisher = {National Academy of Sciences}, + issn = {0027-8424}, + journal = {Proceedings of the National Academy of Sciences} } - @phdthesis{grasinger2019multiscale, - title={Multiscale Modeling and Theoretical Design of Dielectric Elastomers}, - author={Grasinger, Matthew}, - year={2019}, - school={Carnegie Mellon University} + @article{grasinger2021statistical, + title={Statistical mechanics of a dielectric polymer chain in the force ensemble}, + author={Grasinger, Matthew and Dayal, Kaushik and deBotton, Gal and Purohit, Prashant K}, + journal={Journal of the Mechanics and Physics of Solids}, + pages={104658}, + year={2021}, + publisher={Elsevier} } diff --git a/aggregate_mcmc.jl b/scripts/aggregate_mcmc.jl similarity index 100% rename from aggregate_mcmc.jl rename to scripts/aggregate_mcmc.jl diff --git a/scripts/compare_mcmc.jl b/scripts/compare_mcmc.jl new file mode 100644 index 0000000..5542e27 --- /dev/null +++ b/scripts/compare_mcmc.jl @@ -0,0 +1,40 @@ +using Glob; +using DelimitedFiles; +using Quadmath; +using DecFP; + +if length(ARGS) != 3 + println("usage: julia statistical_analysis_mcmc.jl "); + exit(1); +end + +outfile = open(ARGS[1], "w"); +indir = ARGS[2]; +pattern = Glob.GlobMatch(ARGS[3]); + +data_series = Dict(); +for infile in readdir(pattern, indir) + for line in readlines(infile) + var, val = map(strip, split(line, "=")); + val = eval(Meta.parse(val)); + if haskey(data_series, var) + push!(data_series[var], val); + else + data_series[var] = [val]; + end + end +end + +avgs = Dict(); +vars = Dict(); +std_norms = Dict(); +for (k, v) in data_series + avgs[k] = sum(v) / length(v); + vars[k] = sum(map(x -> (x - avgs[k]).^2, v)) / length(v); + std_norms[k] = map(sqrt, vars[k]) ./ avgs[k]; +end +for k in keys(avgs) + @info "$k: μ = $(avgs[k]); σ^2 = $(vars[k]); σ / μ = $(std_norms[k])"; +end + +close(outfile); diff --git a/create_graph_mean_field.jl b/scripts/create_graph_mean_field.jl similarity index 100% rename from create_graph_mean_field.jl rename to scripts/create_graph_mean_field.jl diff --git a/install_dependencies.jl b/scripts/install_dependencies.jl similarity index 100% rename from install_dependencies.jl rename to scripts/install_dependencies.jl diff --git a/scripts/statistical_analysis_mcmc.jl b/scripts/statistical_analysis_mcmc.jl new file mode 100644 index 0000000..5542e27 --- /dev/null +++ b/scripts/statistical_analysis_mcmc.jl @@ -0,0 +1,40 @@ +using Glob; +using DelimitedFiles; +using Quadmath; +using DecFP; + +if length(ARGS) != 3 + println("usage: julia statistical_analysis_mcmc.jl "); + exit(1); +end + +outfile = open(ARGS[1], "w"); +indir = ARGS[2]; +pattern = Glob.GlobMatch(ARGS[3]); + +data_series = Dict(); +for infile in readdir(pattern, indir) + for line in readlines(infile) + var, val = map(strip, split(line, "=")); + val = eval(Meta.parse(val)); + if haskey(data_series, var) + push!(data_series[var], val); + else + data_series[var] = [val]; + end + end +end + +avgs = Dict(); +vars = Dict(); +std_norms = Dict(); +for (k, v) in data_series + avgs[k] = sum(v) / length(v); + vars[k] = sum(map(x -> (x - avgs[k]).^2, v)) / length(v); + std_norms[k] = map(sqrt, vars[k]) ./ avgs[k]; +end +for k in keys(avgs) + @info "$k: μ = $(avgs[k]); σ^2 = $(vars[k]); σ / μ = $(std_norms[k])"; +end + +close(outfile); From 1a100befff19f21b3b3dadfa6a726706066e5860 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Wed, 27 Oct 2021 14:56:25 -0400 Subject: [PATCH 11/17] working on script for post processing --- scripts/compare_mcmc.jl | 75 ++++++++++++++++++++++++++--------------- 1 file changed, 47 insertions(+), 28 deletions(-) diff --git a/scripts/compare_mcmc.jl b/scripts/compare_mcmc.jl index 5542e27..da8c0eb 100644 --- a/scripts/compare_mcmc.jl +++ b/scripts/compare_mcmc.jl @@ -1,40 +1,59 @@ using Glob; using DelimitedFiles; -using Quadmath; -using DecFP; +using GLM; +using DataFrames; +using Logging; -if length(ARGS) != 3 - println("usage: julia statistical_analysis_mcmc.jl "); +if length(ARGS) != 2 + println("usage: julia statistical_analysis_mcmc.jl "); exit(1); end -outfile = open(ARGS[1], "w"); -indir = ARGS[2]; -pattern = Glob.GlobMatch(ARGS[3]); +global_logger(ConsoleLogger(stderr, Logging.Info)); -data_series = Dict(); -for infile in readdir(pattern, indir) - for line in readlines(infile) - var, val = map(strip, split(line, "=")); - val = eval(Meta.parse(val)); - if haskey(data_series, var) - push!(data_series[var], val); - else - data_series[var] = [val]; +indir = ARGS[1]; +subpattern = ARGS[2]; +pattern = Glob.GlobMatch(subpattern * ".out"); + +chunks = map(piece -> split(piece, "-"), split(subpattern, "_")); +param_map = Dict(); +for chunk in chunks + if chunk[1] == "run"; continue; end + param_map[chunk[1]] = Meta.parse(chunk[end]) / 1000.0; +end +@show param_map; + +data_serieses = Dict(); +for subdir in ["standard"; "clustering"; "umbrella"] + data_series = Dict(); + for infile in readdir(pattern, joinpath(indir, subdir)) + for line in readlines(infile) + var, val = map(strip, split(line, "=")); + val = eval(Meta.parse(val)); + if haskey(data_series, var) + push!(data_series[var], val); + else + data_series[var] = [val]; + end end end + data_serieses[subdir] = data_series; end -avgs = Dict(); -vars = Dict(); -std_norms = Dict(); -for (k, v) in data_series - avgs[k] = sum(v) / length(v); - vars[k] = sum(map(x -> (x - avgs[k]).^2, v)) / length(v); - std_norms[k] = map(sqrt, vars[k]) ./ avgs[k]; -end -for k in keys(avgs) - @info "$k: μ = $(avgs[k]); σ^2 = $(vars[k]); σ / μ = $(std_norms[k])"; +Ustar = (param_map["E0"])^2 * abs(param_map["K2"] - param_map["K1"]); +rstar = param_map["b"]; +μstar = (param_map["E0"])*(2*param_map["K2"] + param_map["K1"]); +for (name, data_series) in data_serieses + @show name; + avgs = Dict(); + vars = Dict(); + std_norms = Dict(); + for (k, v) in data_series + avgs[k] = sum(v) / length(v); + vars[k] = sum(map(x -> (x - avgs[k]).^2, v)) / length(v); + std_norms[k] = map(sqrt, vars[k]) ./ avgs[k]; + end + for k in keys(avgs) + @info "$k: μ = $(avgs[k]); σ^2 = $(vars[k]); σ / μ = $(std_norms[k])"; + end end - -close(outfile); From cd793a53195ecfbb34264b348ba0e1e9729e5c50 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Thu, 28 Oct 2021 22:20:04 -0400 Subject: [PATCH 12/17] working version of comparison script --- scripts/compare_mcmc.jl | 42 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 39 insertions(+), 3 deletions(-) diff --git a/scripts/compare_mcmc.jl b/scripts/compare_mcmc.jl index da8c0eb..9b6e3c1 100644 --- a/scripts/compare_mcmc.jl +++ b/scripts/compare_mcmc.jl @@ -3,6 +3,7 @@ using DelimitedFiles; using GLM; using DataFrames; using Logging; +using Printf; if length(ARGS) != 2 println("usage: julia statistical_analysis_mcmc.jl "); @@ -11,6 +12,8 @@ end global_logger(ConsoleLogger(stderr, Logging.Info)); +fmt(x) = @sprintf("%+0.4f", x); + indir = ARGS[1]; subpattern = ARGS[2]; pattern = Glob.GlobMatch(subpattern * ".out"); @@ -18,8 +21,11 @@ pattern = Glob.GlobMatch(subpattern * ".out"); chunks = map(piece -> split(piece, "-"), split(subpattern, "_")); param_map = Dict(); for chunk in chunks - if chunk[1] == "run"; continue; end - param_map[chunk[1]] = Meta.parse(chunk[end]) / 1000.0; + try + param_map[chunk[1]] = Meta.parse(chunk[end]) / 1000.0; + catch e + println("Could not process $(chunk)"); + end end @show param_map; @@ -39,10 +45,15 @@ for subdir in ["standard"; "clustering"; "umbrella"] end data_serieses[subdir] = data_series; end +@show map(keys, values(data_serieses)); +colidxs = Dict(); +for (idx, colname) in enumerate(["step","r1","r2","r3","r1sq","r2sq","r3sq","rsq","p1","p2","p3","p1sq","p2sq","p3sq","psq","U","Usq"]) + colidxs[colname] = idx; +end Ustar = (param_map["E0"])^2 * abs(param_map["K2"] - param_map["K1"]); rstar = param_map["b"]; -μstar = (param_map["E0"])*(2*param_map["K2"] + param_map["K1"]); +pstar = (param_map["E0"])*(2*param_map["K2"] + param_map["K1"]); for (name, data_series) in data_serieses @show name; avgs = Dict(); @@ -56,4 +67,29 @@ for (name, data_series) in data_serieses for k in keys(avgs) @info "$k: μ = $(avgs[k]); σ^2 = $(vars[k]); σ / μ = $(std_norms[k])"; end + rolling_data = map(datafile -> readdlm(datafile, ','; skipstart=1), + readdir(Glob.GlobMatch(subpattern * "_rolling.csv"), + joinpath(indir, name)) + ); + L1s = Dict(); + αs = []; + ϵs = []; + for (varavg, cname, star) in zip([avgs[""][1], avgs[""][3], avgs["

"][1], avgs["

"][3], avgs[""]], + ["r1", "r3", "p1", "p3", "U"], + [rstar, rstar, pstar, pstar, Ustar]) + L1s[cname] = map(i -> begin; + sum(map(rd -> abs(varavg - rd[i, colidxs[cname]]), rolling_data)) / (star*length(rolling_data)); + end, + 1:size(rolling_data[1], 1)); + df = DataFrame(X = map(log, rolling_data[1][:, 1]), + Y = map(log, L1s[cname])); + fit = lm(@formula(Y ~ X), df); + α = coef(fit)[2]; + ϵ = L1s[cname][end]; + println("$cname: $(fmt(α)), $(fmt(ϵ))"); + push!(αs, α); + push!(ϵs, ϵ); + end + println("$(fmt(minimum(αs))); $(fmt(maximum(αs)))"); + println("$(fmt(minimum(ϵs))); $(fmt(maximum(ϵs)))"); end From eb360fba107fe17c7781c09bc5a2341e41a5d9e1 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Fri, 29 Oct 2021 08:30:03 -0400 Subject: [PATCH 13/17] adjusted flipping function to leave phi invariant --- inc/eap_chain.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/inc/eap_chain.jl b/inc/eap_chain.jl index 4320b8e..9130507 100644 --- a/inc/eap_chain.jl +++ b/inc/eap_chain.jl @@ -213,12 +213,17 @@ function flip_n!(chain::EAPChain, idx::Int) move!(chain, idx, π, π - 2*chain.θs[idx]) end +function refl_n!(chain::EAPChain, idx::Int) + move!(chain, idx, 0, π - 2*chain.θs[idx]) +end + pflip_linear(ip::Real) = (1 + ip) / 2; function cluster_flip!(chain::EAPChain, idx::Int; pflip::Function = pflip_linear, ϵflip::Real = 1.0, - ηerr::Real = 1e-10 + ηerr::Real = 1e-10, + flip_f!::Function = refl_n! ) # grow the cluster to the right @@ -258,7 +263,7 @@ function cluster_flip!(chain::EAPChain, idx::Int; if rand() <= ϵflip prev_n̂s = copy(chain.n̂s[:, lower_idx:upper_idx]); for i in lower_idx:upper_idx - flip_n!(chain, i); + flip_f!(chain, i); end if norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) > ηerr From 21479eda6df3b9ea6ea4bddae192cee806b6814c Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Fri, 29 Oct 2021 08:31:20 -0400 Subject: [PATCH 14/17] added more parameter studies --- 2D/run/Ising_2021-10-28.jl | 50 +++++++++++++++++++ ...interacting-dielectric-study_2021-10-27.jl | 50 +++++++++++++++++++ ...interacting-dielectric-study_2021-10-29.jl | 49 ++++++++++++++++++ 3 files changed, 149 insertions(+) create mode 100644 2D/run/Ising_2021-10-28.jl create mode 100644 run/interacting-dielectric-study_2021-10-27.jl create mode 100644 run/interacting-dielectric-study_2021-10-29.jl diff --git a/2D/run/Ising_2021-10-28.jl b/2D/run/Ising_2021-10-28.jl new file mode 100644 index 0000000..1a23616 --- /dev/null +++ b/2D/run/Ising_2021-10-28.jl @@ -0,0 +1,50 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [1e-1; 0.5; 1.0]; +K1s = [0.0]; +K2s = [1.0]; +kTs = vcat(0.25:0.25:1.0, 2.0:1.0:5.0, 10.0:10.0:100.0); +Fzs = [1.0; 5.0; 10.0; 15.0; 20.0; 25.0]; +Fxs = [0.0]; +ns = Int[100]; +bs = [0.5; 1.0; 2.0]; +#bs = [1.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, K1 in K1s, K2 in K2s + push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2, + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b, :run => run)); +end + +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); + +pmap(case -> begin; + + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); diff --git a/run/interacting-dielectric-study_2021-10-27.jl b/run/interacting-dielectric-study_2021-10-27.jl new file mode 100644 index 0000000..659ebd4 --- /dev/null +++ b/run/interacting-dielectric-study_2021-10-27.jl @@ -0,0 +1,50 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [0.0; 1.0; 2.0; 3.0; 5.0]; +K1s = [1.0]; +K2s = [0.0]; +kTs = [1.0]; +Fzs = vcat(0.0:0.05:1.0, 1.5:0.5:5.0); +Fxs = [0.0]; +ns = Int[100]; +#bs = [0.5; 1.0; 2.0]; +bs = [1.0]; +for b in bs, n in ns, Fx in Fxs, Fz in Fzs, kT in kTs, E0 in E0s, K1 in K1s, K2 in K2s + push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2, + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b, :run => run)); +end + +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); + +pmap(case -> begin; + + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); diff --git a/run/interacting-dielectric-study_2021-10-29.jl b/run/interacting-dielectric-study_2021-10-29.jl new file mode 100644 index 0000000..9702a12 --- /dev/null +++ b/run/interacting-dielectric-study_2021-10-29.jl @@ -0,0 +1,49 @@ +println("Hello world!"); + +@everywhere using Printf; + +@show workdir = if length(ARGS) > 0 + ARGS[1]; +else + "."; +end + +@everywhere fmt(x) = @sprintf("%07d", round(Int, 1e3*x)); +@everywhere fmt_int(x) = @sprintf("%03d", x); + +@everywhere function prefix(case) + "E0-$(fmt(case[:E0]))_K1-$(fmt(case[:K1]))_K2-$(fmt(case[:K2]))_kT-$(fmt(case[:kT]))_Fz-$(fmt(case[:Fz]))_Fx-$(fmt(case[:Fx]))_n-$(fmt(case[:n]))_b-$(fmt(case[:b]))"; +end + +cases = Any[]; +E0s = [0.0; 1.0; 3.0; 5.0]; +K1s = [1.0]; +K2s = [0.0]; +kTs = vcat(0.05:0.05:1.0, 1.5:0.5:5.0, 10.0, 100.0); +Fs = [(0.0, 0.0), (0.0, 1.0), (1.0, 0.0), (1.0, 1.0)]; +ns = Int[100]; +#bs = [0.5; 1.0; 2.0]; +bs = [1.0]; +for b in bs, n in ns, (Fx, Fz) in Fs, kT in kTs, E0 in E0s, K1 in K1s, K2 in K2s + push!(cases, Dict(:E0 => E0, :K1 => K1, :K2 => K2, + :kT => kT, :Fz => Fz, + :Fx => Fx, :n => n, :b => b, :run => run)); +end + +@info "total number of cases to run: $(length(cases))"; + +mkpath(workdir); + +pmap(case -> begin; + + outfile = joinpath(workdir, "$(prefix(case)).out"); + if !isfile(outfile) + println("Running case: $case."); + command = `julia -O 3 mcmc_clustering_eap_chain.jl --chain-type dielectric --energy-type Ising -b $(case[:b]) --E0 $(case[:E0]) --K1 $(case[:K1]) --K2 $(case[:K2]) --kT $(case[:kT]) --Fz $(case[:Fz]) --Fx $(case[:Fx]) -n $(case[:n]) --num-steps 1000000 -v 2 --prefix $(joinpath(workdir, prefix(case)))`; + output = read(command, String); + write(outfile, output); + else + println("Case: $case has already been run."); + end + +end, cases); From 633781225f2849d0687177cfa8bb88be1fdc5a62 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Fri, 29 Oct 2021 08:55:35 -0400 Subject: [PATCH 15/17] removed spurious assertion --- inc/eap_chain.jl | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/inc/eap_chain.jl b/inc/eap_chain.jl index 9130507..ce41c84 100644 --- a/inc/eap_chain.jl +++ b/inc/eap_chain.jl @@ -266,17 +266,6 @@ function cluster_flip!(chain::EAPChain, idx::Int; flip_f!(chain, i); end - if norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) > ηerr - println("prev"); - display(prev_n̂s) - println(); - println("curr"); - display(chain.n̂s[:, lower_idx:upper_idx]) - println(); - @show(norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf)) - end - @assert(norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) < ηerr); - new_upper_p = if upper_idx < n(chain) pflip_linear(dot(n̂j(chain, upper_idx), n̂j(chain, upper_idx+1))); else From b61d8f967f7f6fb93662cf64f3aec08f2a3797f6 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Sun, 31 Oct 2021 21:56:44 -0400 Subject: [PATCH 16/17] updated param study --- 2D/run/Ising_2021-10-28.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/2D/run/Ising_2021-10-28.jl b/2D/run/Ising_2021-10-28.jl index 1a23616..e762cf1 100644 --- a/2D/run/Ising_2021-10-28.jl +++ b/2D/run/Ising_2021-10-28.jl @@ -20,8 +20,8 @@ E0s = [1e-1; 0.5; 1.0]; K1s = [0.0]; K2s = [1.0]; kTs = vcat(0.25:0.25:1.0, 2.0:1.0:5.0, 10.0:10.0:100.0); -Fzs = [1.0; 5.0; 10.0; 15.0; 20.0; 25.0]; -Fxs = [0.0]; +Fxs = [1.0; 5.0; 10.0; 15.0; 20.0; 25.0]; +Fzs = [0.0]; ns = Int[100]; bs = [0.5; 1.0; 2.0]; #bs = [1.0]; From 95062a33c44d97be1029692ff6af929dde50fe79 Mon Sep 17 00:00:00 2001 From: Matthew Grasinger Date: Fri, 17 Dec 2021 15:06:07 -0500 Subject: [PATCH 17/17] minor edits to post processing scripts --- scripts/aggregate_mcmc.jl | 31 ++++++++++++++++++++++++++----- scripts/compare_mcmc.jl | 2 +- 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/scripts/aggregate_mcmc.jl b/scripts/aggregate_mcmc.jl index 90f0d94..52bdb1e 100644 --- a/scripts/aggregate_mcmc.jl +++ b/scripts/aggregate_mcmc.jl @@ -3,21 +3,42 @@ using DelimitedFiles; using Quadmath; using DecFP; -if length(ARGS) != 4 - println("usage: julia aggregate.jl "); +if length(ARGS) < 4 + println("usage: julia aggregate.jl [<3D|2D>]"); exit(1); end +dims = if length(ARGS) == 5 + if ARGS[5] == "2D" + 2 + elseif ARGS[5] == "3D" + 3 + else + println("usage: julia aggregate.jl [<3D|2D>]"); + exit(1); + end + else + 3 + end + outfile = open(ARGS[1], "w"); indir = ARGS[2]; pattern = Glob.GlobMatch(ARGS[3]); if ARGS[4] == "dielectric" - writedlm(outfile, ["E0" "K1" "K2" "kT" "Fz" "Fx" "n" "b" "r1" "r2" "r3" "lambda1" "lambda2" "lambda3" "r1sq" "r2sq" "r3sq" "rsquared" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psquared" "U" "Usquared" "AR"], ','); + if dims == 3 + writedlm(outfile, ["E0" "K1" "K2" "kT" "Fz" "Fx" "n" "b" "r1" "r2" "r3" "lambda1" "lambda2" "lambda3" "r1sq" "r2sq" "r3sq" "rsquared" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psquared" "U" "Usquared" "AR"], ','); + else + writedlm(outfile, ["E0" "K1" "K2" "kT" "Fz" "Fx" "n" "b" "r1" "r2" "lambda1" "lambda2" "r1sq" "r2sq" "rsquared" "p1" "p2" "p1sq" "p2sq" "psquared" "U" "Usquared" "AR"], ','); + end elseif ARGS[4] == "polar" - writedlm(outfile, ["E0" "mu" "kT" "Fz" "Fx" "n" "b" "r1" "r2" "r3" "lambda1" "lambda2" "lambda3" "r1sq" "r2sq" "r3sq" "rsquared" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psquared" "U" "Usquared" "AR"], ','); + if dims == 3 + writedlm(outfile, ["E0" "mu" "kT" "Fz" "Fx" "n" "b" "r1" "r2" "r3" "lambda1" "lambda2" "lambda3" "r1sq" "r2sq" "r3sq" "rsquared" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psquared" "U" "Usquared" "AR"], ','); + else + writedlm(outfile, ["E0" "mu" "kT" "Fz" "Fx" "n" "b" "r1" "r2" "lambda1" "lambda2" "r1sq" "r2sq" "rsquared" "p1" "p2" "p1sq" "p2sq" "psquared" "U" "Usquared" "AR"], ','); + end else - println("I don't understand the last argument"); + println("I don't understand the second to last argument"); exit(1); end diff --git a/scripts/compare_mcmc.jl b/scripts/compare_mcmc.jl index 9b6e3c1..c96ed79 100644 --- a/scripts/compare_mcmc.jl +++ b/scripts/compare_mcmc.jl @@ -6,7 +6,7 @@ using Logging; using Printf; if length(ARGS) != 2 - println("usage: julia statistical_analysis_mcmc.jl "); + println("usage: julia compare_mcmc.jl "); exit(1); end