// Native neural diffusion with terminal ensemble-moment calibration.
// Run in a fresh directory; this short run demonstrates mechanics, not research acceptance.
// Fields are nondimensional, with a bimodal stochastic mode and prescribed mass.
import aner.gen
import aner.nn
import aner.tensor
import aner.math
import aner.parallel

nn.configure("cpu", "Float64", "preserve")
let count = 320
let values = DynamicArray<Float64>(capacity: count * 25)
let conditions = DynamicArray<Float64>(capacity: count)
var record = 0
while record < count {
    let condition = (Float64(record % 5) - 2.0) * 0.35
    conditions.append(condition)
    var sign = -1.0
    if parallel.random(1201, record) > 0.5 { sign = 1.0 }
    let fluctuation = 0.4 * (0.9 * sign + math.sqrt(0.19) * parallel.normal(1202, record))
    var row = 0
    while row < 5 {
        var col = 0
        while col < 5 {
            var value = 0.0
            if row > 0 && row < 4 && col > 0 && col < 4 {
                let x = 3.141592653589793 * Float64(col) / 4.0
                let y = 3.141592653589793 * Float64(row) / 4.0
                value = condition * math.sin(x) * math.sin(y)
                    + fluctuation * math.sin(2.0 * x) * math.sin(y)
            }
            values.append(value)
            col = col + 1
        }
        row = row + 1
    }
    record = record + 1
}
let covariates = conditions.to_tensor(rows: count, cols: 1)
let training = gen.field_data(tensor.reshape(values.to_tensor(rows: count, cols: 25), List(count,5,5)),
    covariates, role: "training")
let dx = 0.25
let mass_scale = dx * dx * (1.0 + math.sqrt(2.0)) * (1.0 + math.sqrt(2.0))
let base = gen.diffusion_model(5,5,1,32,64,0.00001,0.2,1211)
let schedule = base.with_diffusion_schedule_features()
let averaged = schedule.with_diffusion_ema(decay:0.999)
let velocity = averaged.with_diffusion_velocity_prediction(loss_weighting: "epsilon")
let spatial = velocity.with_diffusion_sobolev_loss(dx: dx,dy: dx,laplacian_weight: 1.0)
let model = spatial.with_mass_constraint(0,mass_scale,0.0,dx,dx)
let fitted = model.train(training,steps:2000,batch:32,rate:0.001,seed:1213)
// Independent probes repeat the complete training-condition design four times.
// No Gaussian closure: sampling still evaluates the frozen neural denoiser at every step.
let joint = fitted.calibrate_diffusion(training,samples_per_timestep:640,seed:1215,ridge:0.001)
// Keep the learned physical support while regularizing the finite-data amplitude correction.
let calibrated = joint.calibrate_diffusion_neural(training,draws:4*count,seed:1217,strength:0.25)
print(calibrated.loss(training)==fitted.loss(training))
calibrated.save("neural-calibrated.anersci")
let restored = gen.load_scientific("neural-calibrated.anersci")
print(restored.fingerprint()==calibrated.fingerprint())
print(restored.summary())
let inputs = [[0.175]]
let generated = restored.sample_fields(inputs,seed:1223)
print(nn.item(nn.mse(generated,calibrated.sample_fields(inputs,seed:1223))))
print(nn.item(nn.mse(gen.field_mass(generated,dx,dx),nn.scale(inputs,mass_scale))))
print(gen.field_statistics(generated,dx,dx))
// Calibration changes sampled first/second moments. Distribution shape must already
// be learned by the neural chain and independently validated against held-out simulations.
let continued = restored.resume(training,steps:1)
print(continued.steps()==fitted.steps()+1)